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ABSTRACT 

We propose to integrate the Vlasov-Poisson equations giving the evolution of a dynam- 
ical system in phase-space using a continuous set of local basis functions. In practice, 
qq , the method decomposes the density in phase-space into small smooth units having 

fN| ' compact support. We call these small units "clouds" and choose them to be Gaussians 

of elliptical support. Fortunately, the evolution of these clouds in the local potential 
has an analytical solution, that can be used to evolve the whole system during a sig- 
nificant fraction of dynamical time. In the process, the clouds, initially round, change 
shape and get elongated. At some point, the system needs to be remapped on round 
\Q , clouds once again. This remapping can be performed optimally using a small num- 

^sO ' ber of Lucy iterations. The remapped solution can be evolved again with the cloud 

^P , method, and the process can be iterated a large number of times without showing 

significant diffusion. Our numerical experiments show that it is possible to follow the 

2 dimensional phase space distribution during a large number of dynamical times 

i-jH ' with excellent accuracy. The main limitation to this accuracy is the finite size of the 

clouds, which results in coarse graining the structures smaller than the clouds and 



induces small aliasing effects at these scales. However, it is shown in this paper that 
this method is consistent with an adaptive refinement algorithm which allows one to 
~3 ' track the evolution of the finer structure in phase space. It is also shown that the 

generalization of the cloud method to the 4 dimensional and the 6 dimensional phase 
space is quite natural. 
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1 INTRODUCTION 

The problem we are interested in is solving numerically the Vlasov-Poisson system, which writes, in the proper choice of units 



^ + v.V x f-V x <f>.V v f = 0, (1) 



A<f> = 2 f(x,v,t)dv. (2) 



Due to the high dimensionality of phase-space, 2D, where D is the dimension of space, this problem is usually approached with 
the traditional TV-body method, i.e. by approximating the distribution function by a discrete set of particles. However, with 
modern supercomputers, it now becomes possible to start envisaging direct phase-space approach with D = 2 and D — 3. 
In this paper, we are thus interested in solving the Vlasov-Poisson equations directly in phase-space. We consider a new 
implementation in 1 dimension, D = 1, but we shall discuss its extension to higher number of dimensions. We now review 
phase-space methods already studied in the past. After that, we give a sketch of our "clouds" implementation and explain 
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what is new compared to these works. At the end of this introduction, we shall detail the plan of our paper, which is mostly 
devoted to the actual technical details involved in the implementation of our method. 

A fundamental property of Vlasov equation is the Liouville theorem, which states that the phase-space distribution 
function is conserved along trajectories of matter elements in phase space: 

f[x(t),v(t),t] — constant. (3) 

The first numerical methods used in astrophysics to solve Vlasov-Poisson equations in phase space exploited directly this 
property, using the so-called water-bag model (De Packh 1962; Hohl & Feix 1967). The idea of the water-bag model is the 
following. If one assumes that the distribution function is constant within a patch in phase-space, it is enough to follow 
dynamically the boundary of the patch. Numerical implementation of the water-bag model is therefore rather straightforward 
(Roberts & Berk 1967; Hohl & Feix 1967; Janin 1971; Cuperman, Harten & Lecar 1971). Even though this isocontour method 
is quite efficient and accurate, it is in fact very costly. Indeed, the distribution function develops increasing filamentary details 
during evolution by effects of rolling up in phase-space due to differential orbital speeds. Therefore, it is in principle necessary 
to add more and more points to sample the boundary of the patches as time passes. This is one of the major weaknesses of 
the water bag method, which is fine grained by essence, except for initial conditions. 

Other approaches for solving the Vlasov-Poisson equation are grid based, and a large part of the technical developments 
come from plasma physics. One of the most famous numerical implementation, since it inspired many subsequent works, is 
the splitting algorithm of Cheng & Knorr (1976). The splitting scheme consists in exploiting the Liouville theorem in two 
steps, while evolving the system during a time step At: 

f*(x,v) = f{x-vAt/2,v,t), 
f**(x,v) = f*{x,v + V x cj>At), 
f(x,v,t + At) = f**(x-vAt/2,v). (4) 

In the method of Cheng & Knorr, the distribution function is interpolated on a grid using either Fourier method or/and splines. 
It is semi-Lagrangian in the sense that, to compute the value of the distribution function at a grid site, test particles trajectories 
are resolved backwards up to previous time step, where the interpolation is performed. If the original implementation of Cheng 
& Knorr is one dimensional, the generalization to higher number of dimensions is straightforward (e.g. Gagne & Shoucri 1977; 
see Sonnendriicker et al. 1999 for a more recent perspective). The method of Cheng & Knorr was applied first in astrophysics 
by Fujiwara (1981), Watanabe et al. (1981) and Nishida et al. (1981). 

In principle the algorithm of Cheng & Knorr can be used as it is, even when the filamentation effects discussed above 
occur at resolution scale, although it has to be adapted e.g. by using appropriate interpolation procedure to warrant positivity 
of the distribution function and mass conservation (e.g. Besse & Sonnendriicker 2003 for state of the art latest developments). 
However, an elegant solution was proposed by Klimas (1987) to overcome the problem of filamentation. It consists in writing 
the exact equation of evolution of the coarse grained distribution function in velocity space. For that, he assumes a Gaussian 
smoothing window. The modified Vlasov equation giving the evolution of the smoothed distribution function includes a new 
source term. This method was applied to a splitting algorithm using Fourier decomposition (Klimas & Farrell 1994). 

Other grid based methods include hydrodynamic advection schemes: the Lax-Wendroff integration method (Shlosman, 
Hoffman & Shaviv, 1979) or other finite difference methods, using, for the interpolation of the fluxes, standard "upwind" 
and TVD algorithms, convenient to deal with filamentation, such the Van-Leer limited scheme and the piecewise parabolic 
method (PPM), or other schemes such as the flux corrected transport, the flux balanced method (see Arber & Vann 2002 for a 
review and a comparison between these 4 last methods), or the positive and flux conservative method (Filbet, Sonnendriicker 
& Bertrand, 2001). A finite element method was proposed as well (Zaki, Gardner & Boyd 1988), but no further development 
in that direction was performed afterwards, probably because this method involves the inversion of coarse but large matrices, 
a very costly operation in 6 dimensional phase-space. 

It is now worth mentioning two interesting special cases, the Lattice method (e.g. Syer & Tremaine 1995), and the solver of 
Rasio et al. (1989). In the lattice method, the motion of elements of phase space density is restricted to a set of discrete points: 
the time, positions and velocities of "particles" are restricted to integer values, and forces are rounded to nearest integer. 
Such an algorithm has the advantage of solving simplectic equations of motion, which warrants conservation of geometrical 
properties of the system and in the limit of infinite resolution, converges to the "smooth" solution and enforces naturally 
Liouville theorem. The second solver is the spherical code of Rasio et al. (1989), which works in the fully generalistic case. 
The principle of this code is to take full advantage of the perfect knowledge of initial conditions: whenever f(x, v, t) has to 
be determined accurately at some point of space, e.g. to compute accurately the potential, a test particle is followed back in 
time to find its initial position and the value of / associated to it. The needed sampling at each time step is estimated by a 
self-adaptive quadrature routine. This code is therefore quite costly, since at each time step a full set of backward trajectories 
has to be recomputed. It has however the advantage of following with very good accuracy all the details of the distribution 
function and it is probably the most accurate code of this kind available. It presents theoretically the same advantages as the 
water-bag method, but at variance with the later, it is able to preserve smoothness of the distribution function. 
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Finally, alternate ways of solving Vlasov-Poisson equations consist in computing the moments of the phase-space density 
with respect to velocity space and position space, and write partial differential equations for these moments up to some order, 
with some recipe to close the hierarchy (e.g. White 1986; Channell 1995). 

A good numerical implementation of Vlasov-Poisson equations should stick as close as possible to Eq. J3J. In particular, 
it should, as a direct consequence of Liouville theorem, preserve as much as possible the topology of the distribution function. 
More specifically, to render the algorithm TVD (total variation diminishing), and therefore stable, 

• the critical point population i.e. the number of maxima, minima and saddle points of various kinds should be preserved, 
and if it is not possible, should not increase with time. Since all Eulerian implementations use splitting in x and v, this 
condition reduces in ID-ID to preservation of monotonicity in that case. 

• the height of critical points should be preserved, or at worse, the height of local maxima should decrease and the height 
of local minima should increase. This in particular warrants the positivity of the distribution function. 

These conditions state that if the solution deviates from the true one, it should only become smoother. This is the essence 
of modern advection methods, which try to preserve the features of the distribution without adding spurious small scale 
features, such as oscillations around regions with very high gradients. For instance, Van-Leer and PPM methods are TVD, 
at variance with most semi-Lagrangian methods. The higher is the order of the scheme, the more accurate is the resolution 
of fine features of the distribution. However, to preserve the TVD nature of the system, the implementations are lower order 
in sharp transitions regions [condition (a)] and in general nearby local extrema [condition (b)]. This as a result can introduce 
significant diffusive effects. Clearly, we see that to this respect, grid based methods are thus inferior to the water-bag method, 
which optimally fulfills conditions (a) and (b), since it follows isocontours of the phase-space distribution function, and thus 
its topology, by essence. However, the number of sampling points increases with time in the water-bag model and a fair 
comparison with grid based methods should allow adaptive mesh refinement. Clearly, the very particular implementation of 
Rasio et al. (1989) is very successful in that sense. 

The method we propose here is very different from all the works discussed above. The basic idea is to decompose the 
density in phase space into small local units which can be represented using a continuous function with compact support. The 
total density in phase space is the sum of the local functions, to have a fully analytical and continuous representation of the 
phase space density of the system. We will call these small units "clouds" , which will be chosen here to be truncated Gaussians. 
The evolution of the distribution function will be followed by solving the Vlasov equation for each of the individual cloud, 
plunged in the global gravitational potential. We allow the clouds to change shape during runtime, i.e. to transform from 
functions with initial round support (with the appropriate choice of units) into functions with elliptic support. These clouds 
also move: their center of mass position in phase-space follows standard Lagrangian equations of motion, as test particles 
would. If the potential is locally quadratic, our elliptical shape approximation is exact, and the Vlasov equation can be solved 
analytically in that case. It means that as long as the clouds are small compared to the curvature radius of the force, or in 
other words, to the scales of variation of the projected density, this approximation, accurate to second order in space, is very 
good. 

At some point, the local quadraticity condition is not verified anymore, typically after a fraction of orbital time, and the 
analytical solution ceases to be a good approximation of the density of the cloud in phase space: the whole system has to 
be remapped on a new basis of round clouds that give a smooth description of the density in phase space. To resample the 
distribution function with a new set of clouds (including initial conditions), we use a Lucy or Van Citter iterative method. 
This, combined with the fact that remaps are not very frequent, eliminates diffusion almost completely, provided that the 
resolution limit of the simulation has not been reached, i.e. as long as filamentation is not a problem. As a time integrator, 
we use a predictor second order corrector with slowly varying time step. Our algorithm is thus fully second order in time and 
in space, and nearly simplectic, since in the case the time step is constant, it reduces to standard Leapfrog. 

As a matter of fact and as we shall see, our method deals very well with filamentation, as it naturally coarse grains the 
distribution function at small scales. It is quite close to finite element methods, except that the elements are of changing shape 
and that the sampling grid moves with time. It is also close in some sense to the water-bag method, but in a coarse grained 
way, as by construction equation 10 is verified exactly as long as the potential is locally quadratic. The remap procedure is 
iterative and costly, that is one of the drawbacks of the method, similarly as in finite difference methods. However, we do not 
need to perform it too often (typically every 10-20 time step). The Lagrangian nature of the method and the reinterpolation 
scheme makes our method very weakly diffusive, but it is not TVD: aliasing effects can appear in regions with high curvature. 
However, for the numerical cases we studied, these effects are not critical. Positivity of the distribution is enforced with 
Lucy reconstruction, but not with Van-Citter. Note that our method can be used only in warm cases, i.e. in cases where the 
distribution function is smooth and has some width in velocity space. It is not appropriate for the cold case, e.g. to describe 
in phase-space formation of large scale structure in standard cosmological models. 

This paper is organized as follows. In § 2, we present the theoretical background intrinsic to our method, i.e. study the 
dynamics of a phase-space cloud in a quadratic gravitational potential. We perform a full perturbative stability study, taking 
into account small deviations from quadraticity. In § 3, we discuss the actual numerical implementation. In § 4, we test our 



4 C. Alard & S. Colombi 

code by performing simulations of a stationary solution, an initially Gaussian profile and an apodized water-bag (a top hat 
in phase-space). In § 5, we propose an adaptive refinement procedure, which allows one to increase resolution where needed, 
if details of the distribution function need to be followed at smaller scales. Finally, § 6 summarizes the results and discusses 
some perspectives for the method, in particular its extension to higher number of dimensions and the treatment of the cold 
case. 



2 THE METHOD: CONCEPTS 

2.1 One dimensional equations 

In this section we will show that for a quadratic potential, the Vlasov equation for a cloud has an analytical solution. Let 
4>{x,t) be the gravitational potential of the system. Locally the quadratic approximation will read: 

<j>(x,t) = ao(t) x + ai(t) x + a 2 (t). (5) 

By taking the following cloud equation it is possible to show that the Vlasov equation forms a closed system: 

f(x,V,t) = G [\o(x,t) + \i(x,t) v + \ 2 (t) v 2 ] , (6) 

where G is any smooth function (with continuous derivatives). The functions Xi(x, t), i = 0, ..., 2, which determine the geometry 
of the cloud in phase space can be obtained by solving the Vlasov equation: 

2l + v 2l _^ ftf =0 . (7) 

dt dx dx dv 
It is interesting at this point to separate the general motion of the cloud from the evolution of its internal geometry. We 
will thus use Lagrangian coordinates, which are defined by (x* = x — xa,v* = V — Vo), where (xq, vg) are the phase space 
coordinates of the center of gravity of the cloud. Assuming that the cloud phase space density f(x, v, i) become f(x*,V*,t) in 
this referential, we can write: 

f + !£,._ W |L = o, (8) 

at ox ov 

fix* ,v*,t) = G [a (x* , t) + Ai (x* ,t)v* + A 2 it) u* 2 ] , (9) 

^2- = -2ao(t)xG - ai(t). (10) 

By inserting Eq. @ in Eq. @, we obtain a quadratic polynomial in v* , which has to cancel for any value of v* . Hence, 
we obtain 3 equations: 

d 

—A ix*,t)-2A 1 ix*,t)aoit)x* = 0, (11) 

— Ai(x*,t) + — A (x*,i)-4A 2 (i) ao(t) x* = 0, (12) 

d A 2 it) + JL;Aiix*,t) = 0. (13) 



dt dx 

The general solution of Eq. i1K'>l) is 



dt 
By inserting the solution for Ai in Eq. 1121 and solving for Ao (£*,£) we obtain: 

,2 dip 



Ai(x*,t) = -^x + Mt)- (14) 



A (a; ,t) = - 



d A 2 . , . 

+ 4A 2 (t)a (t) 



.'• --j^-x* + ip2(t). (15) 



dt 2 

And finally by reporting Eqs. 11 H and 1151 into Eq. 11 U . we obtain a quadratic polynomial in x*, which must be zero for any 
value of a;*. The zeroth order of this polynomial gives ^2 =constant. The first order of this polynomial reads dtpi/dt = —2aoipi. 
However, in the referential of the center of gravity, the first moments of / with respect to x* and v* should be zero. This 
implies that the form Ao(x*,i)+Ai(a;*,i)u*-|-A2(i)w* 2 , that we know now to be a polynomial of order 2 in x* and v* with time 
dependent coefficients, does not have any term either in x* and v*\ it is a quadratic form in x* and v* , implying ipiit) = 0. 
The second order of this polynomial writes 

1 d A2 , , dA 2 „dan , , . .„„. 

Eq. p6^ can be solved with respect to ao(t): 

a (t) A 2 (t) 2 = -i I ^lA 2 it)dt + C. (17) 



A cloudy Vlasov solution. 5 

Since 

f d 3 A 2 , , , , d 2 A 2 , , , 1 fdA 2 \ 2 , N 



we can write 

i>'-i(7)' t '*'* , - <1 <19 > 

It is possible to simplify this equation a little more by using the following substitution: A 2 (£) = w(t) 2 (here we assume that 
A 2 (£), which corresponds to the velocity dispersion, is positive). We obtain an equation which is easier to use for numerical 
integration: 

5^ + »(*)*(«) = ^- ( 2 °) 

We can now finally rewrite Eq. J3J as follows 

f(x*,v*,t) = G\/3o(t) a;* 2 +/3i(i) x^*+/3 2 (t)w* 2 +/3 3 ] , (21) 

with 

/*>(*) = i^+2A 2 (t)ao(t), (22) 

*(*) - -^, (23) 

/3 2 (t) = Aa(t), (24) 

and /?3 can be set to zero without any loss of generality. We then see that Eq. I|190 has a simple geometric interpretation, 
since it can now be rewritten 

0o(t)0 a (t)-~fc(t) a =2C, (25) 

which says that the area of the elliptical section in phase-space defined by Eq. 112 1 1 , proportional to 

A = -kA = = (26) 

remains constant with time. 

2.2 Dynamical properties of the clouds 

In this section, we develop a number of simple analytical calculations which will help to implement the numerical simulations. 
These analytical models will also help us to understand the limits of our approach. 

2.2.1 Small amplitude oscillations 

The evolution of the cloud is given by equation Eq. 1)200 . In general, the solution of this highly non-linear equation can be 
approached only by numerical means. However, if we consider small amplitude oscillations around an equilibrium position, 
it is possible to linearize the equations and find a simple analytical solution. We should consider a stationary equilibrium 
position with a quadratic potential given by: otoit) = qo- In this case, the solution to Eq. 121)1 is: 

w(t)=w =(-\ . (27) 

Note that due to the Poisson equation, 

it is possible relate the stationary potential qo to the mean density po around the cloud: 

90 = Po- (29) 

Now, we introduce a time dependent perturbation of this stationary regime: 

ao(t) — qo + eqi(t) and w(t) — Wq + ewi(t). (30) 

Then Eq. J2SI reads: 

^-^- + ±qoW 1 {t)+w q 1 {t)=Q. (31) 
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A general solution to Eq. 1331 can be obtained by using a Fourier transform. We define the following Fourier transforms: 

wi(w) = / Wi(t) e llJt dt and Qi(w) = / qi{t) e 1 ^ dt. (32) 

V 2n J \l2ix J 



Using these definitions, the general solution of Eq. I|33|l writes: 

Wilu) = 2u>o —5 o + ^0 d( w ~ u o) + fci <5( w + w o)> (33) 

W — o;q 

where fco and k\ are two arbitrary constants and 

w = 2^2^). (34) 

Eq. llMMI shows clearly that a resonance may occur between the potential and the cloud oscillations at the frequency uio- 
However for this resonance to occur effectively, it would be required that all the nearby clouds which contribute to the local 
density resonate also at this frequency. But, the resonant frequency of the nearby clouds can be the same only if the local 
density around these clouds is identical, which would require a constant density whatever the position. A resonance happens 
only in this case, and we will discuss this very special case later. The dynamical properties of the cloud are also dictated by 
the motion of their center of gravity which is given by Eq. (1101 . The solution to Eq. ill I It can be decomposed as the solution 
of the homogeneous equation and a given solution to the whole equation. The homogeneous equation reads: 

^+2a (t)x G (t) = 0. (35) 

Setting XG(t) — xa + exi(t) and noting that Eq. 1351 1 is similar to Eq. 12011 . we find the general solution in the linearized 
regime: 

Xi(w) = 2w —^ — r, + fco S(uj - lu ) + fei 6(u + Wo), (36) 

with a resonant frequency, a>o = wo/2. Thus all the dynamical frequencies of the clouds are multiples of the fundamental 
frequency \/2qo, and thus, with the help of Eq. 1291 . we infer that the dynamical time of the cloud is: 

Tcloud — — =■. (37) 



Obviously this dynamical time has little meaning for a cloud experiencing large potential variations on a short time scale. 
In such case the non-linearity is dominant and the linear approximation is not suitable. However, for most clouds in a given 
system, this approximation will give an estimate of the local time-scale. 



2.2.2 Resonant self- oscillations 

As noted in the former section, resonances between the cloud motion and oscillations would require that the frequency Wo is 
constant. Using Eq. I29|l and Eq. I|34fl . we see that this condition Wo = constant would require that the density is constant. 
The simplest example of a system having constant projected density is a line in phase space, with a density constant along the 
line, the oscillation frequency for any point on the line is the same whatever the position angle of the line in phase space. Thus 
such system can have self-oscillation, and out of this case there is no regime of self-oscillation of the cloud system. Note that 
this one-dimensional system is similar to the case of the sphere of constant density that has been shown to have self-pulsation 
in 3 dimensions. 



2.2.3 Deviation from the quadratic approximation 

The cloud approximation assumes that the potential at the scale of the cloud is quadratic. Even if the quadratic approximation 
is a good description of the local potential, when evolving the system, we expect that the errors to the quadratic approximation 
will accumulate, and that at some point the cloud model will deviate from the proper solution. Our numerical scheme should 
stop to evolve the system before the deviation is large, and this is the purpose of this section to estimate the magnitude of 
this deviation as a function of the dynamical time describing the system evolution. We will consider an additional cubic term 
in the potential, and we will study the effect of this term in the linear perturbative regime. The new equation for the potential 
reads: 

4>(x*,t) =a (t) x* 2 +e/3(t) x" 3 . (38) 

We introduce the corresponding perturbation for the density in phase space: 

f(x*,v*,t) = f (x*,v*,t) + ef 1 (x*,v*,t). (39) 
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Using the Vlasov equation [Eq. 0], we obtain the following equations: 



dfo df * , 9/o , . 

-TTT + tt— V - 2 ao X ■ H — = 0, (40) 

at ox* ov* 

dfi.dfi , . dfi * 2 9/o n ,. n 

^ + ^7- "2«o^ ^-3/3. — =0. (41) 

We already know the solution to Eq. 1401 from section 2.1. Since clouds having gaussian profile are particularly interesting 
for numerical applications, we make the following choice in Eq. JUJ 

fo(x*,v*,t) =exp[Ao(a;*,i)+A 1 (x*,t) v* + A 2 (t) V* 2 ] . (42) 

Now, it is possible to find a general solution to Eq. 1411 by using the following functional for /1: 

/1 (x*,V*,t) = /o (x*,w*,t) X 770 (t) + rfi (t)x* 3 + 772 {t)x* 2 v* +rj 3 (t)x*v* 2 + 774 (t)v* 3 \ . (43) 

By combining Eq. 1431 with Eq. (I4H . we obtain a polynomial in v* of order 3. Since the coefficients of this polynomial must 
be identically zero, we find the following system of differential equations: 

% + 277 2 -6a»74 = 0, (45) 

at 

%--6/3A 2 + 377i-4Qr?3 = 0, (46) 

at 

^--2ar ?2 -3/?A 1 =0. (47) 

at 

2.2.4 Effect of time dependent non-quadratic terms: a practical case 

In practice the potential felt by a cloud moving in a given system can be quite different from quadratic. The Poisson equation, 
Eq. 12811 . shows that the non-quadratic terms are related to local gradients in the density distribution. The accurate estimation 
of the local potential depends on the details of the density distribution, thus to quantify the effects we will have to adopt a 
given density distribution. However, in general, as the cloud moves, the system evolves and changes its potential, but this 
potential variation is slow and can be neglected for our practical purpose. Let us assume that the density distribution is given 
by the family of profiles: 



P(x) = Q (f ) , (48) 



where h is a scale-height. 

Using the properties of this density distribution and the Poisson equation, it is easy to analyze the scale properties of the 
cloud equation. The coefficients of the local potential ao and (3 in Eq. I|38|l can be rewritten using the scaled variable x* : 

1 d d> , „, , *. , . 

a " " 2 !>x*~ =p{x ^ = q (* )' ( 49 ) 

1 dp{x*) = l_ dq(x*) 

3 dx* Sh dx* 

Eq. I|49|l shows that ao is independent of the scale of the density distribution. As a consequence, Eq. I|19|l is also scale 
independent, and thus the quadratic solution A 2 (t) is unaffected by the scale of the density distribution. But the behavior 
of the coefficient f3 is different. According to Eq. 15UH . ft is scale dependent. It is possible to evaluate the effect of this scale 
dependence by forming an equation for 774 (£) similar to Eq. 1191 for Aa(f). By combining Eqs. 1441 . 145H . 1461 . and 147H . it is 
possible to construct the following equation for 774(4): 

a d0 . dA 2 d 2 a da dn 4 d 2 7? 4 1 d 4 rj 4 % . , 

- 6 d i - A2 - 15/3 ^r + 3 ^^ 7?4 + 10 ^r^F + 10Q0 ^ + 2^ + 18a ° 7?4 = a (51) 

It is easy to notice that the scaling property of (3 and the scale invariance of «o and A2 imply that 774 scales like 1/h. Which 
means that the correction introduced by the cubic term is inversely proportional to scale. Thus, we may solve the equation 
for 774 at a given scale and generalize the result to any scale using the former scaling rule. The solution for a fixed scale can 
be calculated using Eqs. I|49fl . 1500 and 151|l . 

For our numerical implementation, we will adopt a fixed gaussian density distribution for the system, q(x) = exp(— x 2 ). 
Given initial conditions in phase space, the evolution of the cloud crossing this density distribution can be solved numerically 
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Figure 1. The maximum deviation from the quadratic approximation as a function of the initial position xq. The distribution has a 
scale length h ten times larger than the cloud. The coordinate xg has been normalized by the scale length of the distribution h. Note 
that the maximum deviation is at about 3h from the center of the distribution. 
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Figure 2. The shape of the deviation due to the cubic term as a function of the strength of the cubic term. We chose an initial position 
xo which maximizes the deviation from the quadratic approximation (xo ~ 3h). Note that in case the distribution has a smaller scale 
length the error must be rescalcd accordingly. In practice, the minimum size of the distribution is about 2 times the cloud scale length, 
for which the deviation can be 5 times larger. 



by integrating our former system of equations. The motion of the center of gravity of the cloud itself can be computed easily 
by estimating the motion of a point mass in the potential of the system. The orbit of the cloud is defined by its initial position 
xo, and its initial velocity. We will study the case of zero initial velocity and variable position xo. Initially the cloud will be a 
round gaussian of velocity dispersion equal to 1. The numerical integration of these equations shows that the deviation from 
the quadratic approximation is maximal for an initial position close to 3h (sec Fig. 1). The maximum deviation observed in 
the system is directly related to the requirement to perform a Lagrangian remap. Thus we should study the behavior of the 
error to the quadratic approximation near the 3h initial position. The relevant plot is shown in Fig. 2. As discussed before, 
the curve showing the deviation due to the cubic term is generic. The deviation for another scale length h of the density 
distribution can be obtained by rescaling this curve. Since the scaling is inversely proportional to h, the error to the quadratic 
approximation will be dominated by the crossing of the shorter scale density structures. 
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3 THE METHOD: ALGORITHM 



In this section, we study practical implementations of the method. We examined two approaches, one assuming constant 
resolution in phase-space, that will correspond to what we call a "cloud in mesh" (CM) code, the other one allowing local 
refinement in phase-space using adaptive refinement trees, that we therefore call tree-code. There are several issues to be 
addressed while implementing the method. We list them here in the same order as they will be treated below: 

• Phase space sampling of the distribution function ("S I3.1I I : the question here is to decide how to chose our set of clouds such 
that it reproduces at best a given distribution function in phase-space. This is necessary not only to set up initial conditions, 
but also to resample the distribution function during run time with a new set of round clouds. Indeed, we know from our 
perturbative analysis (§ 12.20 that deviations from local quadratic behavior of the potential increase with time, and that at 
some point the clouds will have wrong axis ratio (be too elongated) and wrong orientation. To (re)sample the distribution 
function, we propose to use Gaussian clouds located on a (possibly adaptive) grid, with their masses estimated using either 
Van-Citter or Lucy deconvolution algorithms. 

• Solving the Poisson equation ("§ 13.21 : the issue here is to to estimate accurately the forces exerted on each cloud as well 
as errors on their determination. These latter will indeed be used to quantify deviations from a local quadratic behavior of 
the potential. This is fairly easy in our one dimensional problem, since the force exerted on a point of space is simply given 
by the difference between the mass at the right of the point and the mass at the left of the point. However, we will try here to 
experiment methods which can be, in principle, easily generalized to higher number of dimension without any serious increase 
in complexity. 

• Run time implementation and diagnostics (§ 13.31 : the choice of the time step implementation is important. Here we 
propose a predictor corrector of second order, which makes our method full second order both in space and in time. We 
discuss how to compute the slowly varying time step associated to this integrator. We also examine when a remap of the 
distribution with a new set of round clouds has to be performed, by using a criterion based on the cumulated errors on the 
forces due to non local quadratic behavior of the potential. 

3.1 Phase space sampling of the distribution function: practical implementation 

When starting from given initial conditions, f(x,v), it is necessary to set up an ensemble of round clouds that all-together 
reproduce at best /. During runtime, these clouds get more and more elongated and the approximation for a local quadratic 
potential breaks: a remapping of the distribution function is needed with a new set of round clouds. The measured distribution 
function at this time becomes new initial conditions that, again, have to be sampled accurately. Here we explain in detail how 
we proceed to perform this remapping. In § 13. 1.11 we discuss the cloud shape and the mean cloud inter-spacing. Our choice 
is a Gaussian cloud of typical radius R, truncated at 47? and separated from its nearest neighbors by a distance \/2R. In 
8 13. 1.21 we explain how we compute the cloud masses in order to reproduce at best the distribution function. The methods 
proposed are Van-Citter and Lucy deconvolutions algorithms. In § 13.1.31 we propose a weighting scheme for filtering small 
scale noise in the estimate of the distribution function obtained from the clouds. Finally § 13.1.41 examines possibilities of 
enforcing conservation of basic quantities such as total momentum and total mass without adding significant diffusion. In 
addition to that, Appendix[X]details the algorithms used to compute the distribution function from the clouds. It shows how 
to speed up the calculation in the case the clouds are round, thanks to the separability of the Gaussian window. 

3.1.1 Choice of cloud shape and spacing 

To sample the distribution function in phase space at a given time, we use an ensemble of clouds which are initially round 
(this can be always true if an appropriate normalization for positions and velocities is set). In this section, we assume that 
they are disposed on a rectangular grid of spacing A g . The choice of our cloud shape is entirely determined by specifying 
function G [Eqs. @ and Q21|l ]. which should present the following features: 

(i) It should have compact support. Indeed, we need the cloud to be the least extended as possible, since we use a local 
quadratic expansion of the potential within the cloud to compute the forces exerted on it. 

(ii) It should be sufficiently smooth in order to resample a distribution function with continuous first and second derivatives 
(e.g. for refinement procedure as discussed later). This is an essential feature of our method and a condition for it to perform 
well. 

A good choice for function G is a truncated exponential, which has the great advantage of being separable and makes the 
cloud Gaussian. Smoothness condition (ii) forces us to truncate this cloud rather far away from its center. Our practical 
choice is a 4 sigma cut-off, -R max = 4.R. 1 This still induces small discontinuities of the order of 3.10 -4 , as illustrated by right 

1 More exactly i? ma x = 3.95-R, to avoid the circle of radius i? ma x intersecting any of the cloud center. 
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panel of Fig. [3] To minimize their effect we will use a weighted estimator to compute /, as detailed in § 13.1.31 The remote 
nature of the cut-off has another disadvantage, which is one of the main drawbacks of our method, when it will be extended 
to higher number of dimensions: a large number of clouds contribute to the sampled distribution function at a given point of 
phase-space, as illustrated by lower right panel of Fig. 

Finally, we have to determine the radius R (corresponding to a one sigma deviation) of our cloud as a function of 
grid spacing for best sampling the distribution function. Basically, the choice of R determines by how much the resampled 
distribution function, /, will deviate from the true /. In order to sample smoothly the variations of / over the sampling grid, 
R/A g should be of the order of unity. To minimize the cost of the sampling, it should be kept as small as possible. To find 
R, we estimate the quadratic error a due to the representation of space by a finite number of Gaussian functions. Since our 
approximation requires / to be nearly flat at the scale of the cloud, we will evaluate the residual in the case where the function 
to represent is constant, / = 1. For an infinite flat distribution a reads 



S-oo [f( x > v )- f( x ,v] 



2 

dxdv 



f dxdv 

J — oo 

For gaussian clouds, / is represented by the functional: 

_[(x-iA g ) 2 + (v-jA 
2-kR? 



A 2 



./] 



2R 2 



Note that the representation f(x, v) can be rewritten as a convolution using Dirac series 



(52) 



(53) 



^' u ) = rfr£/ ex P 



2-kR? 

i,3 



(uf + ul) 



2R? 



8(ui + x — i A g )5(u2 + v — j A g ) du\du2. (54) 



This rewriting as a convolution suggests that the equations should be analyzed using Fourier transform. Since both the 
numerator and denominator in Eq. 1521 are the norm of a function in real space, the transformation to Fourier space can be 
done easily by using Parseval theorem. Using the symbol / to represent the Fourier transform of /, Parseval theorem reads: 

f°° r ~ i 2 f°° I - = I 2 

/ \f(x,v) — f(x,v)\ dxdv = \f(k x ,k v )~f(k x ,k v )\ dk x dk v . (55) 

J — oo J — oo 

And with the help of the convolution theorem: 

/(^,fc„)=^exp[-2 7r 2 i? 2 (k 2 x + k 2 v )]S(k x -i/A g )S{k v ^j/A g ). (56) 

ij 

Thanks to the rapid fall of the exponential, this sum is well approximated by a small number of terms. In particular a good 
approximation at the typical scale of interest is to consider only the terms (i,j) verifying \i\ + \j\ < 1: 

f \f(k x ,k v )- f(k x ,k v )\ 2 dk x dk v ~4exp[-(2nR/A g ) 2 ] S(0) 2 . (57) 

The divergent factor 5(0) disappears once we apply the normalization by the denominator in Eq. 1521 . We thus finally find: 

a ~ 4 exp [- (2-KR/Ag) 2 ] . (58) 

Our choice R — v2/2 corresponds to a — 10 -4 , consistent with the discontinuities brought by our 4 sigma cut-off for the 
gaussian. This is well illustrated by lower-left panel of Fig. |31 which displays the deviations of / from unity for our choice 
of R/A g . Experimentally, we measure that their mean square is indeed equal to 10 -4 , while the difference between local 
minima to local maxima is four times larger. Notice again that the truncature at 7? ma x adds a source of noise at very small 
scales which has the bad property of presenting a significant skewness in its distribution. We shall see in § 13.1.31 that with 
appropriate reinterpolation of function /, it is possible to remove these sources of noise to a great extent, in fact totally in 
the case / =constant considered here. 

3.1.2 Choice of deconvolution scheme 

Once the cloud shape is chosen as well as the size of the sampling grid, the problem of finding the cloud masses in order to 
sample correctly the distribution function remains open. To do that, an iterative procedure is necessary. We adopted in this 
paper two simple algorithms, Van-Citter and Lucy as we discuss now. 
The sampled distribution function can be written 

f(x,v) = ±^2MiG(x-xh,v-vh), (59) 

i 

where Mi is the mass of the cloud and function G(x,v) is normalized in such a way that its integral is V . The quantity V is 
given by V — A X A„, where A x and A„ correspond to grid site inter-spacing distance along x and v coordinate respectively. 
It is the area associated to each cloud so that f f(x, v)dxdv — J^\ Mi. Note that we used implicitly the following notations: 

G(x,y) = G(p t) x 2 + f3 1 xy + f3 2 y 2 ). (60) 
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Figure 3. Sampling of f(x,v) = constant = 1 by truncated Gaussians. We consider a set of 5 X 5 clouds in a periodic square of size 
L = 5. With these units, each Gaussian is of size R = V2/2 and is truncated at i? ma x = 3.95H. Upper left panel shows the variations 
of the sampled distribution function / if the clouds would have infinite extension (darker and lighter regions correspond to / > 1 and 
/ < 1 respectively). They range with a good accuracy from 2<r ~ 2.1 X 10~ 4 to — 2<r [Eq. 1581 1. Upper right panel is the same, but the 
truncation is applied, which explains the discontinuities observed in /. In that case, fluctuations of / are larger and are not symmetrically 
distributed, ranging approximately from —3.8 X 10~ 4 to 2.3 X 10 -4 . Lower right panel shows the cloud counts, i.e. for each point of space 
the number of cloud contributing to it. This count varies significantly, between 21 (light) and 27 (dark). Finally, lower left panel displays 
the values of f(x, v = x) (the diagonal of the images) for the untruncated (thick smooth curve) and truncated case (thin irregular curve). 
The truncation adds a significant source of noise with a significant skewness in its distribution, which can become a source of systematic 
effects. This effect has to be minimized with the appropriate weighting scheme. 



We have to find Mi such that 

f(%G,vh) =fi = /(kg,«g) = fi, 
for each cloud position {xq,Vq) in phase-space. A good first guess is simply 

M,° = fiV. 



(61) 



(62) 



However, with that choice of masses, / will basically be equal to the convolution of / with function G. Some deconvolution 
algorithm has to be applied in order to fit better function /. 
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The Van-Citter algorithm consists simply as follows. Given /" as computed at iteration n, the iteration n + 1 writes, 

o t — }i - Ji , 
M n+1 = M n + gny_ (gg) 

The algorithm is applied until some convergence criterion is fulfilled, max(|5™|) < S c . To maintain the domain of calculation 
as compact as possible during runtime, only values of /; verifying fi > 5 m ln are taken into account. In the simulations 
shown below, we take 5 m ln — S c /2, but <5 m i n should in principle be kept as small as possible to have best conservation of 
the moments of the phase-space distribution function, e.g. energy. Similarly, after convergence, clouds contributing little to 
the reconstruction are depleted: one uses reciprocal of Eq. il(i2t to estimate approximately their contribution, i.e. clouds with 
\Mi/V\ < 5min arc eliminated. This procedure has the defect of augmenting slightly aliasing effects in the neighborhood of 
the regions where / cancels. 



Given the sources of uncertainties associated to our choice of function G, 5 C should be large enough compared to 10 /, 



max • 



where / max is the maximum value of /, since at this point we would capture spurious features and convergence would be 
made difficult due to the discontinuities brought by the cut-off. Experience shows that <5 c // max ranging from 0.0005 to 0.002 
is a good compromise. Also, convergence might become difficult due to other aliasing effects. As well shall see later, function 
/ builds finer and finer structures with time, which cannot be reproduced correctly by our mapping when they appear at 
scales smaller than ~ R. Convergence is rendered very difficult in that case: one might chose therefore to stop iterating when 
n reaches some value, typically of the order of 10 according to our practical experiments. If convergence is not reached at 
this point, it is clear anyway that fine structures of function / wont be reproduced correctly, even with a large number of 
iterations. 

The major defect of Van-Citter algorithm is that it does not warrant positivity of the distribution function. An alternative 
to Van-Citter is Lucy deconvolution algorithm: 

s? = fi/fr-i, 

M? +1 = M?J2(l + 5 3) G ( x 'G-x^vh-v 3 G ). (64) 

3 

Such an algorithm does not only warrant positivity, it is also in principle more accurate than Van-Citter, since the error 8" 
is relative instead of absolute. 2 In practice, Lucy does better than Van-Citter, but is a little more unstable and is also more 
subject to aliasing effects. For our fixed resolution simulations, we adopted Lucy algorithm, but we used Van-Citter when 
testing refinement, as discussed later. 

To illustrate the methods, Fig.|I]shows some results obtained when trying to resample the following distribution function, 
used later as initial conditions for some of our simulations: it is a top-hat apodized with a cosine: 

f{x,v) — p, x + v < TL , 

f(x,v) = ^ {COS [| ( Vz 2 + V 2 - K) /ft apo ] + l} , X 2 +V 2 <ft + 2ft apo , (65) 

where p — 1, TZ = 0.7 and 7J. ap o = 0.3. Function f(x,v) is shown and plotted in top panels of the figure. Middle panels 
of the figure are the same but the difference 8f — f — f between the reconstructed distribution function and the true one 
is considered, when a small number of clouds is used, with inter-spacing A g = A^ = A„ = 0.2. We used 10 iterations in 
Eqs. J(».'~!> and H64fl and a cut-off value 5 m i n = 0.5 x 10~ 4 . A larger number of iterations would not change significantly the 
level of convergence. Since A g is of the same order of 7?. ap o, aliasing effects are rather significant, particularly for Lucy where 
they reach a few percent magnitude, because of the positivity constraint intrinsic to this algorithm. Lower panels of the figure 
are the same as middle panels, but a larger number of clouds was used, with inter-spacing A g = 0.05 now small compared 
to 7£apo- As a result, aliasing effects are much less significant, of the order of 0.001 — 0.002 and the difference between Lucy 
and Van-Citter has decreased. This shows that, to sample correctly variations of / over some length scale I, we need A g to 
be sufficiently small compared to I, typically l/A s ~ a few unities. Reversely, we see that in runtime, the resampled / can be 
fully trusted only at coarse graining scales of the order a few A g 's. 

3.1.3 Filtering small scale noise 

After using either Van-Citter or Lucy algorithms, one obtains from interpolation 1591 a function / which by definition 
reproduces at best the true values of / at the sampling points {x 1 q,Vq), given some convergence criteria. We just discussed 
problems of aliasing, which are intrinsic to the method and cannot be really avoided. They can be only reduced by using 
e.g. refinement procedures discussed later, or just by decreasing cloud inter-spacing. However, it is possible to reduce defects 
discussed in § 13.1.11 namely small scale variations due to the finiteness of phase-space sampling (upper left panel of Fig. [3J 

2 In that case, 5 C is not anymore expressed in term of /max- 
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Figure 4. Comparison of Van-Citter and Lucy algorithms for reconstructing a distribution in phase-space. Upper left panel displays 
the function f(x,v) we aim to reproduce, as given by Eq. Iti5l [Darker regions correspond to larger values of f(x,v)]. Upper right panel 
shows f(x, 0) as a function of x. Middle left panel displays Sf = f — f, where / is obtained by using Van-Citter algorithm. A very similar 
result would be obtained for Lucy. The middle right panel gives Sf(x, 0) as a function of x, both for Van-Citter (thin black curve) and 
Lucy (thick grey curve). To sample the distribution function, an inter-spacing of A g = 0.2 and a cut-off -R max = S.QSA^ was used. The 
lower-left and lower-right panels are the same but for higher resolution, A g = 0.05. Note that, for the four lower panels, the weighting 
scheme discussed in 5 13. 1.31 has been used to estimate /, following Eq. 1661 . 
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Figure 5. The difference Sf(x, 0) as a function of x, between reconstructed (with Van-Citter) and true distribution function, whether 
Eq. 1591 (thin black curve) or Eq. 1661 is used to compute /. The function / to sample is the same as in Fig. Eland the resolution is the 
same as in lower panels of this figure. 



and discontinuities due to the truncature of the clouds (upper right panel of Fig. |3J , to warrant smoothness of the sampled 
distribution function. This is indeed necessary to avoid propagating this intrinsic noise during runtime, when a large number 
of reinterpolations is performed. 

To minimize these defects, we perform a weighting different from Eq. 1591 1 : 



f(x,v) = 



(66) 



where W = J^ fc , G(A x k, A v l) ~ 1 to ensure proper normalization. This weighting procedure is nearly equivalent to Eq. 159H . 
since we expect the function '}2 i G(x — x r G ,v — Vq) to present small variations with (x,v) as shown by upper right panel of 
Fig. However, it has the advantage of smoothing quite efficiently small scale noise as illustrated by Fig. 03 For instance, 
in regions where / =constant, it recovers exactly the true value of / within these regions (with a possible offset due to 
uncertainties in the convergence of the reconstruction). Hence, our weighting scheme would give exactly / = 1 everywhere in 

Fig.H 

In what follows, when a new mapping with round clouds is necessary to resample an existing cloud distribution, we shall 
use Eq. I66H to compute / from the old set of clouds and Eq. 1591 to compute / from the new set of clouds in iterative 
procedures I63H and 164H . We decided to do so in the last case to avoid introducing a bias in the new cloud masses estimate. 



3.1.4 Enforcing basic conservations 

An intrinsic property of deconvolutions methods such as Lucy's or Van-Citter's is that they conserve very well mass. However, 
we will have in practice to perform many remaps of the distribution function: any small but systematic deviation from 
mass conservation might have, on the long term, dramatic consequences. As a result, it might be necessary to enforce mass 
conservation, but this is not easy without adding diffusion. In order to reduce as much as possible this latter we proceed as 
follows. Suppose that the total mass of the system should be equal to M, and let be M = J^ Mi the total mass obtained 
from the reconstructed clouds. If M > M, we compute the cumulated positive mass residues: 



AAf 4 



Y, 5M ^ 



(67) 



i,« last >0 
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where 5 Mi — 8\ aBt V for Van Citter algorithm and 8 Mi = Mi<5| ast for Lucy algorithm: these quantities correspond to the 
remaining uncertainty on the mass determination for each cloud after last iteration, n =last. Then, for each cloud having 
^last y q^ we j ncrease j^ s mass ^ a f ac tor a 5Mi with a = (M — M)/AM+. We proceed similarly if M < M, but by considering 

AM- = J2 dM i- ( 68 ) 

i,<5> a3t <0 

Our mass conservation scheme should minimize as much as possible diffusion effects, since the correction calculated for each 
cloud is proportional to the uncertainty on its mass determination. Furthermore, only clouds for which the mass should be 
increased (or decreased) in the direction of M — M are modified. Note that if a > 1, it means that our correction will tend 
to be above the residues. We can clearly fear for diffusion in that case. Once total mass conservation is taken care of, we can 
pay attention to total momentum, 

P v =Y,Mivh, (69) 

i 

that we always force to be zero by appropriate corrections of the velocities, if needed, since there is no external force exerted 
on the system we consider. 

3.2 Solving the one dimensional Poisson equation. 

In this section, we explain in details how we solve Poisson equation. To do that, one needs first to estimate local projected 
density, which is rather simple with our Gaussian clouds, as discussed in § 13.2.11 Then one has to estimate the force exerted 
on each cloud, as well as its slope. To do that, we propose two methods, a tree-code approach based on decomposition of space 
on a binary tree (§ 13.2,211 and a "cloud in mesh" (CM) approach based on sampling of space with a regular grid f§ 13.2.31 . 

3.2.1 Calculation of projected density 

The projected density can be obtained for any value of x by summing all the individual contributions of the clouds. For a 
given Gaussian cloud, with elliptic basis and arbitrary orientation, integration over velocity space of Eq. 1211 is in fact simple. 
With the help of Eq. 11911 . we obtain 



p ^ t) = M \Rw)^{-^T) x ) (70) 

for a cloud of mass M. In this equation, we thus neglect the truncature of the cloud over velocity space, but this should not 
have any significant consequence. Indeed, in practice, to make sure that the total projected mass of each truncated cloud is 
equal to its true mass, we renormalize Eq. I|7()|l by a factor very close to unity. 

3.2.2 Practical implementation: a tree code. 

Keeping in mind that we want to generalize the code to 2D and 3D, we experimented an implementation which can be, in 
principle, easily extended to higher number of dimensions. The Poisson equation simply writes, in the appropriate units 

The approximation taken in this paper thus consists in simply assuming that p = constant within a cloud. The force, 4 given 
by 

—fT= I P(y,t)dy~ I p(y,t)dy = M ri &t-Mutt, (72) 

is then locally fitted by a straight line within the cloud. 

To estimate the force, we decompose the x axis hierarchically on a binary-tree (Fig.|SJ, by dividing successively segments 
into two parts of equal length, until each projected cloud intersects with at least N s tree cells. Of course, there is an uncertainty 
on the real number of intersecting cells, Mater > JV S . It should be at least N s > 2 and in fact large compared to that if we 
want to estimate errors on the forces. To compute the force within a cloud, given the number iVinter of tree cells intersecting 
with its projection on x axis, we store (i) the force Fi exerted on each cell center, Xc, and (ii) a weight ua proportional 
to p(xc — Xg) — p{ x c % ) as given by Eq. 17011 . This weight is of course more important when the cell center is close to the 

3 Also, it can be possible but rather unlikely that, e.g., M > M and AM+ = 0. It is not worth enforcing mass conservation in that case, 
since the reconstruction itself is probably already very biased due e.g. to strong aliasing effects. 

4 In this paper, we indistinctly liken the quantity — V</> to a force or an acceleration. 
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Figure 6. Sketch of our tree structure used to compute the forces within each cloud. Here, we consider the case N B = 3, i.e. each cloud 
must intersect with at least 3 cells of the tree at a given refinement level. Note that for computing the forces for cloud A, only cells 7,8,9 
and 10 of level 3 of refinement will be used, even if more refinement is available due to cloud B. 

projected cloud center. Then, given a list of (xi,Fi), i = 1, iVinter, we use a simple weighted least square fit to adjust the force 
by a local straight line Fa t (x) = —2ctox — qi [see Eq. JjJ]. Our estimator for the error on the force then reads 



AF 
~FZ~ 



Winter 

— VtOi \Fi-F m (xh)]' 



(73) 



max/ M 2 £f=l tK V~< ,_, 

This weighted quantity roughly quantifies the deviation of the force from a straight line within the cloud. It is renormalized 
by the maximum force, F max , which in the units chosen here is equal to the total mass M of the system. The error on the 
force is indirectly related to the projected size of the cloud along x axis. If the cloud gets considerably elongated, the error 
on the force will become larger in a fixed potential. Therefore, one might use A_F/.F max < some small value as a criterion to 
decide whether a remap of the distribution function with a new set of round clouds is necessary. However, as we shall discuss 
more in detail in § 13.3.31 what really counts is the cumulated error, 

(74) 



£ cum = / AF(t)dt. 



To make sure that the force is estimated consistently the same way for all the clouds (in particular to preserve approximately 
local smoothness properties) , we stop refining the tree locally as soon as A/inter > N e , even if the tree has many branches due 
to much smaller clouds at the same location x, as shown by Fig. |5] Since we chose here _R m ax = 47? for the truncature of 
the cloud, one might think that iV s = 8 should be taken, to have at least two sampling cells per typical length-scale, 7?, of 
the clouds. However, we noticed that 7V S = 2 is enough to have a very accurate determination of the force with our choice of 
R/A x = v2/2. This is illustrated by Fig.[7| To estimate the instantaneous error on the force, a larger value of N B is necessary, 
typically iV s > 8, while N s > 4 is enough to compute the cumulated error as used in § 13.3.31 



3.2.3 Alternate practical implementation: "CM" code 

The results presented in this paper use a tree-code. The tree-code has the advantage of being rather flexible and easily 
generalizable to higher number of dimensions. It also represents a natural ground for adaptive refinement in phase-space, as 
we shall see later. To simplify the approach, it is also possible to assume fixed resolution in space, which allows one to use 
FFT or relaxations methods to solve Poisson equation, similarly as in particle-in-mesh codes (PM). However, at variance 
with PM, a much more accurate "cloud-in-mesh" interpolation (CM) is performed. In our ID case, the resolution of Poisson 
equation is rather trivial. The projected density of the clouds is simply calculated on a grid of inter-spacing A s in such a 
way that the number of intersecting grid sites is always larger than N s for any cloud. The force is then estimated exactly the 
same way as discussed in S 13.2.21 In the current implementation, this fixed resolution code is much faster than the tree-code, 
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Figure 7. The calculated force for the examples considered in Fig. [3] The left and the right panels correspond to the case A g = 0.2 
and A g = 0.05 respectively. On each panel, there is a smooth thick grey curve: it corresponds to the exact force derived directly from 
integrating Eq. 1651 . There is a thin step like dotted curve: it corresponds to the force computed for each cell of the binary tree at 
the resolution needed to sample each cloud with a least JV S cells. Here we took N s = 2, and in practice we obtain -/V; ntor = 2 or 3. 
There are stars nearly perfectly fitting the smooth thick grey curve: they give the force calculated at the center of each selected cloud 
for doing the figure (a subset of all the clouds). A segment passes through each of these stars: it gives our weighted least square fit 
F&t{ x ) = —2aox — ai, which trivially reduces to the line passing through the two sampling points when -/Vj ntor = 2. It is interesting to 
notice that the estimate we obtain for the force at the cloud centers is very accurate, although one might notice a slight offset between 
symbols and the thick smooth grey line on left panel. The projected length of each segment on x axis is equal to the projected length 
of each cloud (up to cut-off). As we see, significant deviations from our local quadratic assumption are expected for the potential on 
left panel. These deviations are much less visible on the right panel. The typical weighted deviation calculated from Eq. 1731 by using 
N B = 8 gives AF/F ~ 0.03 and 0.002 for left and right panel respectively. 

especially at the moment of remap, where a fast convolution method can be used to estimate the phase space distribution 
function (Appendix IA2L However, note that the tree-code part can still be considerably optimized for the remap part, since 
the fast convolution method can in principle be generalized to a non structured grid. 

3.3 Time stepping implementation and diagnostics 

In this section, we discuss our run time implementation: a second order predictor corrector, as described in § 13.3.11 This 
makes our approach second order both in time and in space. The determination of the slowly varying time step is discussed 
in § 13.3.21 as well as other diagnostics, such as energy conservation. Finally, § 13.3.31 examines the critical issue of deciding 
when performing a new sampling of the distribution function with a set of round clouds, by studying cumulated errors on the 
determination of the forces. 



3.3.1 Time stepping implementation: global predictor corrector 

In our approach, we need a number of parameters to describe completely a set of clouds i, i — 1, . . . , N to t, interacting through 
gravity. These parameters can be chosen as follows: 

- the mass of the cloud, Mc this one does not change with time, except when a remap with a new set of round clouds is 
performed; 

- the position of the center of the cloud, x x G (t); 

- the velocity of the center of the cloud, VQ(t); 

- the acceleration of the center of the cloud, that we write Fi(t); 

- the parameter cto(t) in Eq. (|^J, which is nothing but the projected density at position [(EgC0>' 1 'g(*)]j 

- the function Wi(t) = \/Aj(t); 

- the time derivative of function Wi(t): Dv)i(i) = dwi/dt; 
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- the area of the elliptical section of the cloud divided by n: Ai — Ai/n: this one does not change with time, except when a 
remap with a new set of round clouds is performed. 

Recall that the parameters /3q, P\ and f3 2 defining the shape of the cloud are entirely determined by function u>i, its time 
derivative Dwi, a l and Ai, namely, from Eqs. 12211 . 123H and 1241 . with the help of Eq. 1201 : 

/3g = {Dw t f + — - pi=-2wiDwi, f3 2 = (w 1 ) 2 . (75) 

[AiWiy 

To evolve the clouds we use standard second-order predictor corrector algorithm, which is known to preserve quite well 
simplecticity, a feature essential in phase-space. To simplify the algorithm, we take the same global time step for all the 
clouds. Our run time implementation can be split up into six main parts: 

(a) Predictor step: 

:eg(Wi/2) = xh{tn) + -dtnvh{tn), (76) 

Wi(t n+1/2 ) = w l (t n ) + -dt n Dwi{t n ) . (77) 

(b) Force calculation: the force is calculated as explained in § 13.2.21 at the "predicted" positions x l G (t n+1 f 2 ) within each 
cloud of "predicted" projected density given by Eq. 170t using A2(i n +i/2) = {u>i(t n+ i/ 2 )] 2 ■ This gives, for each cloud, the 
parameters Fi(t n+1/2 ) and cto(t n+1/2 ). 

(c) Corrector step: 

xh(t n +i) = x' G (t n+1/2 ) + -dt n v G (t n ) + -dt 2 l F i (t n+1/2 ), (78) 

vh(tn+i) = v Q (t n ) +dt„Fi(t n+1/2 ), (79) 



1 1 2 

Wi(t n+1 ) = Wi(t n+1 / 2 ) + -dt n Dwi(t n ) + -dt n 



A 2 wf(t n+1/2 ) 



2ao(t n+1 / 2 )wi(t n+1 / 2 ) 



Duii(t n+1 ) = Dwi(t n ) + dt n 



2a Q (t n+1 / 2 )wi(t n+1 / 2 ) 



(80) 
(81) 



_Alwt(i n+1/2 ) 

(d) Outputs: various quantities such as all the informations on the clouds, the density in phase space, etc, can be output 
at this point. 

(e) Diagnostics: calculation of next time step, test for Lagrange remap, energy conservation check: as we are going to 
describe more in detail below, we use a time-step constrained by the slope of the force, —2a % {t n+ i/ 2 ). The fact that the next 
time step, dt„+i, depends on a quantity calculated half a time step before does not affect significantly the quasi-simplecticity 
of our integrator. The important thing here is that dtn should vary slowly with time. At this point, we also test if it is 
necessary to remap the distribution function with a new set of round clouds, due to the accumulation of deviations from local 
harmonicity of the potential, as studied in § 13.3.31 

(f ) Lagrange remap, if needed, as explained in § 13.11 

Note importantly, thus, that our algorithm is second order both in space and in time. It is quasi-simplectic in the sense that 
it would be time reversible, equivalent to generalized leap frog, if the time step was constant. 

3.3.2 Diagnostics: "Courant" condition and energy conservation 

There is a list of diagnostics to perform during runtime. The most important are those related to time stepping. What matters 
in our Lagrangian approach is that orbits in phase-space should be sampled with sufficient number of points as well as cloud 
shape variations during the trajectory. From the analysis of S 12.2.11 that means, using Eq. I|37|l . 

dt < — = , (82) 

maxi \f 2a l 

where C is a "Courant" parameter small compared to unity. In practice we find that we should have C<0.01. For all the 
simulations made in this paper, we took C = 0.01. 

Another diagnostic is of course total energy conservation, which writes, in our units, 



1 [ 2 1 

£ tot = - / v f(x,v,t) dx dv - - 



d 4X~M 2 

dx 



dx = constant. (83) 



This equation assumes that the system is finite of total mass M. We added a term proportional to M to compensate for the 
divergence in the integral J(d(f>/dx) dx so that the constant in Eq. I|83|l is finite. 

One might use instantaneous deviations from energy conservation to constrain the time step, but we chose here to test 
energy conservation only as a consistency check. We shall indeed see in §|l]that energy is very well conserved in our code. 
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Finally notice that mass is conserved by definition, except when a resampling of the distribution is performed with a new 
set of round clouds. During this step, mass conservation is enforced as explained in § 13. 1.41 as well as total momentum. 

3.3.3 Lagrange remap 

The assumption for a quadratic local potential is in general always violated to some extent during runtime since it corresponds 
to variations in the projected density p(x) f8 l2.2.4H . What counts, as already mentioned in S 13.2.21 is the cumulated error on 
the force, which translates in clouds with increasingly wrong ellipticity parameters. This effect, even if always small during a 
time-step (see 8 12.2.40 . cumulates with time. In particular, the clouds get inevitably spuriously elongated and some "hairy" 
structure appears during evolution (Fig.[2j. This is why, we have at some point to remap the distribution function with a new 
set of round clouds, an operation that we call "Lagrange remap". 

To illustrate this point, we consider here the same example as in 8 12.2.41 the evolution of a distribution function initially 
Gaussian, i.e. given by 

ft \ - ( lx 2 +V 2 \ 2 . 2 ^ — 2 

f{x,V) = pexpl-- — — — J, x +v <TZ , 

f(x,v) = ^p exp (~ x + V J {cos [| (\Ac 2 +v 2 - Tl) /fcapo] + l} , x 2 + v 2 <ft + 27l ap o. (84) 

This initial profile will be used for the set of simulations studied in § ^ The same apodization than in Eq. I|65|l is performed 
to regularize the function at its edges: Eq. 1841 approaches Gaussian initial conditions only if h <C 1Z. For the simulation 
considered here, we take h — 0.2, p — 1/(2tv1i 2 ), 1Z—1, 7£ a po = 0.1. The cloud size is R — h/10 — 0.02, similarly as in 8 12.2.41 
and the inter-cloud spacing is thus given by A g = \/2R — 0.028. To compute the force and the error on it [Eq. (I73H 1. we take 
N s = 8 as advocated in 8 13.2.21 

Left panel of Fig.|H]gives the maximum instantaneous error on the force as a function of initial position, x ln i. It shows that 
the deviations from local quadratic approximation globally increase with time, as expected. What is important here, though, 
is the accumulation, E cum , of little kicks at each time-step due these deviations, as given by Eq. 1741 . This quantity, shown 
as a function of time in middle panel of Fig. |H| is homogeneous to a velocity. It should remain small compared to velocity 
resolution: 

pmax 

=£^ < E max « 1. (85) 

If condition 1851 is violated, a new sampling of the distribution function with round clouds should be performed. We notice 
that a typical value of the threshold, -E max = 0.05, i.e. a 5 percent cumulated error, corresponds, as expected, to a fraction 
of dynamical time, t ~ 0.4. Interestingly, but not surprisingly, if condition 1851 is enforced, the number n r ema P of time-steps 
between each remap is pretty stable, as shown on right panel of Fig. |H| For our Courant parameter choice C — 0.01 and for 
£" max = 0.05 we find n rcmap ~ 15. Its variations decrease with time, as expected, while the system converges to a stationary 
state. As a result, it is a very good approximation to keep n romap fixed, which is quite useful, since it is not necessary in that 
case to estimate the cumulated error, and as a consequence, N s — 2 can be taken to speed up considerably the calculation of 
the force, without any significant loss in accuracy. 

One would thus be tempted to fix n re map as a function of our Courant condition using the following formula, 

n rcmap (C) = j^-y (86) 

with, typically, no.oi — 15 to produce a cumulated error on the force of the order of 5 percent. This estimate is only valid 
for the initial conditions considered here, but according to arguments of 8 12.2.41 the same result should be roughly found for 
other configurations: with C ~ 0.01, a Lagrange remap should be performed every 10-15 time-steps. 

One might be worried that a cumulated error of a few percent per Lagrange remap is too much: accumulation of errors 
after a few thousand remaps might still be too large. However, on has to take into account the fact that the cumulated error, 
as given by Eq. 1741 . does not take into account possible global cancellations during the dynamics. In particular, if a given 
symmetric potential is fixed, we see that a cloud plunged in the right side of the potential will be the object of next to 
quadratic distortions opposite to those in the other side: they will cancel each other, which makes the solution rather stable 
on the long term, even if our absolute cumulated error becomes very large. This is illustrated by Fig. which compares the 
results obtained for the simulation of Fig. |H1 without remapping at t — 6 to the exact solution. The "hairy" effect is quite 
visible, but the solution is rather close to the correct answer at the coarse level. 



4 SOME APPLICATIONS 

In this section, we show how our tree-code performs for three different kinds of initial conditions. Similar results should be 
obtained for the CM code, so we do not feel necessary to show them here. 
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Figure 8. Left panel: the maximum instantaneous error, AF/F miK , on force calculation [Eq. 1731 1 as a function of initial cloud position, 
measured in a simulation without remapping. The initial conditions considered here are a truncated Gaussian of size "R. as given by 
Eq. 1841 . The initial position, x ln j, of the cloud is expressed in units of h. To reproduce a case similar to the one studied in Fig. the 
cloud size R is given by R/h = 0.1. Three dynamical times are considered, t = 0.5 (solid), 1.5 (dots) and 6.0 (dashes). It is difficult to 
make an exact comparison to Fig. given that fact that our estimate of AF/F max is a quadratic dispersion. However, the general order 
of magnitude found on the error agrees with the results of § 12.2.41 The dotted curve which corresponds approximately, dynamically, 
to the parameters chosen to make Fig. agrees well, at least qualitatively with Fig. Middle panel: the maximum cumulated error, 
E™/A,, as a function of time. A remapping should be performed each time this error exceeds some threshold. Right panel: the number 
of time-steps, n rem ap, between each remap as a function of remap count, when the condition 1851 is fulfilled, with £ max = 0.05. The 
simulation used to make this plot is exactly the same as the one used for left and middle panels, except that Lagrange remap was 
enforced. The simulation was stopped after 3000 time steps, at t ~ 39, corresponding to 35 dynamical times according to Eq. 1371 with 
po = P- Note that the variations of n rem ap decrease with time: n re map seems to converge progressively to some fixed value, as expected, 
while the system relaxes to a stationary state. 



The first model, considered in § 14.11 is a stationary profile, verifying df/dt = 0. This is a crucial test: the code should 
be able to maintain a stationary profile during many dynamical times to a very good accuracy. We shall see that our code 
passes this test with great success. 

The second model, considered in § 14.21 is a Gaussian distribution function: this kind of very smooth initial conditions 
is well adapted to our method, that should be able to follow their evolution quite accurately until details appearing during 
runtime become too small to be resolved by our sampling clouds. We shall see that at this point our numerical solution still 
gives a very good coarse level version of the full resolution one. 

The third model, considered in § 14.31 is a top hat apodized by a cosine. The evolution of such a kind of distribution 
function allows us to test directly to which extent our code can maintain a region with / ^constant. In particular, effects of 
aliasing can be quantified easily. 

For each of the models considered above, we perform one "low resolution" simulation with A x — A v = 0.02, and one "high 
resolution" simulation with A x = A„ = 0.005. We run these simulations using a Courant condition C — 0.01 and a Lagrange 
remap every 15 time steps, following the conclusions of the analyses above. To perform the resampling, we use Lucy algorithm 
with 10 iterations exactly, whatever the residue obtained between the resampled distribution and the previous one, following 
discussions of § 13.1.21 To maintain the domain of contributing values of / finite, we proceed as explained in the paragraph 
just after Eq. 1631 with <5 m i n = 5.10 -5 . In all the cases, we aim to evolve the system during approximately 30 dynamical times 
(corresponding to 60 orbital times) , a goal that we achieve in the following way: we set up initial conditions such that po is of 
order two in Eq. 1371 for the top-hat and the Gaussian case, and evolve them up to t — 100. For the stationary solution, the 
choice of po is of order unity and we evolve the simulation up to t — 150. With our Courant condition, this amounts finally to 
approximately 7000 time steps for the stationary simulation and 8000 time steps for the simulations with Gaussian and top 
hat initial conditions, corresponding to a large number of remaps, about 470 and 530, respectively. 

This way we are able to quantify the limitations of our approach, namely coarse graining and aliasing effects due to the 
finite size of the sampling clouds at the moment of remap. We shall see that, as expected, the low-resolution simulations give 
a very accurate coarse grained version of the high resolution ones. 
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Figure 9. Apparition of a hairy structure after a few dynamical times. The simulation considered is the same as the one as in the two 
left panels of Fig.[S] at t = 6. The top left panel gives the "correct" answer, i.e. f(x, v) obtained by remapping the density field with new 
round clouds using condition 1851 with -E max = 0.05. The top right panel shows the "wrong" answer obtained without remapping. The 
bottom panel displays the difference between left and right panels, which shows the hairy structure appearing due to the divergence of 
the axis ratios of the clouds. 



4.1 First application: stationary profile 

A family of stationary solutions is given by (Spitzer, 1942) 

f s (x,v) =p 
with 



cosh I — 



lv" 

" >:i>i -2^ 



2-Kp 



(87) 
(88) 
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Figure 10. Left panel: test of energy conservation for a stationary solution: The quantity (Etot — £tot theo)/£tot theo is displayed as a 
function of time, where Stot is measured in our simulation according to Eq. 1831 and £ tot thco is given by Eq. 1901 . The parameters used 
to do the simulations are p = 2 and a v = 0.2. The thick and the thin curves correspond to the high and the low resolution simulation, 
respectively. Right panel: deviations from stationarity: the quantity max \f(x,v, t) — f B (x,v)\ is represented as a function of time, where / 
is the simulated / and f B (x,v) is given by Eq. 1871 with additional apodization as explained in the text. The continuous and the dotted 
curve correspond to the high and the low resolution simulations, respectively. 



For further reference, the total mass writes 



and the total energy [Eq. 1831 1 reads 



M = 2(2tt 



1/4 P 1/2 *1 /2 , 



(89) 



£tot=K^) 1 P 1 ^'- (90) 

To perform the simulations, we set up initial conditions with same apodization as in Eqs. 1651 and 1841 . to make the support 
of the function compact. To achieve an initial set up sufficiently close to the stationary solution, we must have a v and 
g x sufficiently small compared to apodization radius, 1Z. Our choice here is 72. = 1, 72. a po = 0.2, p — 2, a v — 0.2 [hence, 
Ox — 0.2 ~ <j v ]. With this set up, the number of clouds contributing during runtime is about 80000-120000 and 4200-7500 in 
the high and the low resolution simulation respectively. This number decreases slowly with time as a result of the truncation 
procedure used to maintain the support of / compact (see S 13.1.21 . In principle we could fine tune this procedure in order to 
maintain the number of sampling clouds approximately constant with time, but we decided not to do so at this point, because 
it affects only the tails of the distribution function. 

Fig. 1101 examines energy conservation and deviations from the stationary solution as functions of time. As expected, 
energy conservation is very good both for the high and for the low resolution simulation, better than 0.5 percent. It however 
worsens with time. This is mainly a consequence of our truncation of the phase-space distribution function at <5 m i n = 5.10 -5 . 
Note that the maximum difference between the true distribution function and the simulated one tends to augment linearly 
with time at a certain point, even for the high resolution simulation, as illustrated by the right panel of Fig. I1UI For the 
high resolution simulation, this is again mainly due to the truncature of the tails of /. Indeed, the difference 5f between 
the numerical and the analytical solutions measured at the center of the fluctuation, (x,v) = (0,0), is only of the order of 
— 2.10 -4 for t = 150. For the low resolution simulation, an other effect adds up to the truncation of the tails, namely a weak 
diffusion at the maximum of / at the moment of Lagrange remap. When fluctuations become of the order of the cloud size, 
one can indeed expect two competing effects, according to the details of the shape of the sampled distribution function: (i) 
a coarse graining effect, where small scale fluctuations tend to be progressively smeared out, (ii) an aliasing effect already 
discussed in § 13.1.21 where, on the contrary, artificial contrasts are created with the appearance of spurious oscillations. Since 
the distribution function we sample here is smooth, effect (ii) is not present. However, although a deconvolution method such 
as Lucy's (or Van-Citter's) aims to minimize effect (i), it cannot reduce it completely. This effect is very small, since it induces 
a rather weak diffusion, even for the low resolution simulation, of the order of 2 percent only at the center of the fluctuation 
at t = 150. 

Fig. 1111 shows the zeroth and second moments of f(x, v) with respect to velocity, namely the projected density and 
the velocity dispersion, as functions of x, at latest output, t = 150. Even though the low resolution simulation has already 
significantly diffused, with a two percent depletion on the maximum of /, the overall solution is still pretty close to the 
correct one. The deviations observed in the tails both in the high and the low resolution simulation are a consequence of the 
truncation procedure used to maintain the computing domain finite. 



A cloudy Vlasov solution. 23 



10.0000 




1.0000 


r 1 


>N 


: / \ : 


U) 


/ \ 


5 0.1000 

"D 


: / \ 1 


"a 


■ / \ ■ 


« 0.0100 

Q_ 


; / \ I 


0.0010 


' \ \'\ 


0.0001 





10.0000 



1.0000 




c 
o 



<u 0.1000 

Q_ 



o 
o 

> 



•1.0 -0.5 0.0 0.5 1.0 
x 



1.0 -0.5 0.0 
x 



0.5 1.0 



Figure 11. Mean projected density, p(x) = J f(x,v)dv (left panel) and velocity dispersion, ii rms (i), v 2 Ins (x) = {Sv 2 ) = J f(x,v)v 2 dv — 
[J f(x, v)vdv] 2 (right panel), as functions of position, for the stationary solution examined in Fig. I1UI The solid curve corresponds to 
Eq. 1871 with apodization as described in the text, which explains the bending of the tails visible in the right panel. The dotted and the 
dashed curves give the results obtained from the high and the low resolution simulations, respectively. 



4.2 Second application: a gaussian as initial conditions 

The gaussian initial conditions we consider are specified by Eq. 1841 . with 1Z = 0.8, 7J. ap o = 0.2, p = 4 and h — 0.2. With 
this set up, the number of clouds contributing is roughly 95000-120000 and 6600-7700 for the high and the low resolution 
simulation respectively. 

Fig. 1121 shows the distribution function in phase-space at various times, t — 0, 10, 40 and 100. As expected, a spiral 
structure appears, which roles up with time. At some point, it becomes so thin that it disappears due to finite resolution. Of 
course, this event happens earlier for the low resolution simulation. Visually, this latter seems to represent a very good coarse 
grained version of the high resolution simulation. This can be examined more in details on Fig. 1131 which displays projected 
density and velocity dispersion at the same instants as in Fig. 1121 At t — 10, it is not possible to distinguish yet between low 
and high resolution. The difference between the two simulations is the most significant for t = 40: in that case, details are 
lost in the low resolution simulation, but it gives a rather accurate coarse grained version of the high resolution one. At late 
time, t — 100, details have nearly disappeared even in the high resolution simulation and the agreement between low and high 
resolution remains excellent. 

Fig. 1141 displays the relative deviation from energy conservation as a function of time. For the high resolution simulation, 
results are very similar to what was obtained for the stationary solution: total energy decreases slowly with time, mainly 
because of our truncature of the tails of the phase-space distribution function, but is conserved with a precision better than 
0.5 percent at t = 100. The behavior observed for the low-resolution simulation is more complex. Indeed, energy first increases 
up to t — 50, due to coarse graining effects which introduce a slight bias in the tails of the distribution function. At variance 
with the high-resolution simulation, these effects dominate over those due to the truncature until details in the distribution 
have completely disappeared. Note that the maximum of /, that should stay constant, presents small variations between initial 
and final time. The net result is an increase of 0.1 and 2 percent for the high and the low resolution simulation, respectively. 
This is probably due to aliasing effects, which are however not expected to affect significantly energy conservation. 

Examination of the second part of Fig. 1131 suggests that the system is relaxing towards a stationary regime. Fig. 1151 
compares the phase space distribution function at the latest stage of the high resolution simulation to a stationary solution 
given by Eq. 1871 . To compute the parameters of the stationary solution, we use Eqs. (I89H and 1901 . For the total energy, we 
take the initial value obtained from the (apodized) Gaussian profile and obtain p ~ 3.79 and a v — 0.298, which corresponds 
to the dotted curve. It globally agrees well with the simulated solution (thick grey curve), if ones keeps initial conditions as 
a reference for comparison (solid thin curve). Note that the agreement is improved in the high / part if p is changed to 4 
(dashed curve), but remains unperfect: the dashed curve is not able to reproduce the clear transition between a plateau at 
low values of / and the central belt shape, which suggests a more complex stationary regime, with two components. 
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Figure 12. Phase space distribution function at various times, for Gaussian initial conditions. Left and right panels correspond to the 
high and the low resolution simulation, respectively. Darker regions correspond to higher values of /. 
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Figure 13. Mean projected density (left panels) and velocity dispersion (right panels), as functions of position, for the simulations with 
Gaussian initial conditions. Each row of panels corresponds to a given time. The solid curve and the dotted curve are for the high and 
the low resolution simulation, respectively, except for t = 0: in that case, the curves are calculated according to Eq. 1841 . 



4.3 Third application: a top hat as initial conditions 

To set up initial conditions, we use Eq. 1651 with TL = 0.8, 7J. apo = 0.2 and p — 1. The number of clouds contributing to / 
varies between 190000 and 350000 and between 11000 and 29000 for the high and the low resolution simulation, respectively. 
Figs. Hoi and 1171 show the resulting distribution function in phase space and its zeroth and second moments, at various 
times, t — 0, 10, 40 and 100, both for the low and the high resolution simulation. The systems builds with time a two 
component structure: most of the region f — p remains in a compact structure with a roughly elliptic shape rotating and 
pulsating around the center, while a ring appears around it by an effect of rolling up. The appearance of this ring, formed of 
a very fine spiral structure transporting a small fraction of the total mass (see Fig. I17B . explains the significant increase of the 
number of clouds contributing to the distribution function. At t — 100, our pulsating core+ring structure has not converged 
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Fig . 1131 (continued) . 



to any stationary solution, as expected, but remains topologically stable on the coarse level, if one takes into account the fact 
that mass migrates slowly from the edges of the central patch to the external ring. 

A careful examination of Figs. 1161 and 1171 shows that the agreement between low and high resolution is pretty good, 
although not as impressive as for the Gaussian case: some details differ even at the coarse level, particularly in the region of 
transition between the external ring and the central patch. These differences do not affect significantly energy conservation, 
as illustrated by Fig. 1181 in both simulations, energy decreases slowly with time, mainly as a result of our truncation of the 
tails of the phase-space distribution function, and is conserved with an accuracy better than 0.2 percent. 

Note finally that aliasing effects seem to be more significant here than for the Gaussian simulation, although still quite 
reasonable: at t — 100, we find for both simulations that the maximum of / is about 1-2 percent larger than it should be. 
This happens at the edges of the central patch. At the center, we measure /(0, 0) ~ 1.005 and 1.0007 for the low and the high 
resolution simulation, respectively, which represents a very small deviation from unity. 
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Figure 14. Energy conservation for the simulations with Gaussian initial conditions. The relative deviation from energy conservation is 
displayed as a function of time for the high resolution simulation (thick curve) and the low resolution simulation (thin curve): £tot and 
£tot ini correspond to total energy and initial total energy, respectively, as measured in the simulations according to Eq. 1831 . 
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Figure 15. Comparison of the simulation with Gaussian initial conditions with a stationary solution. Function f(x, 0) is displayed in 
the left panel, while function /(0, v) is displayed in the right panel. On each panel, the simulation (thick grey curve) is compared to the 
stationary solution given by Eq. 1871 with same energy and same mass (dots). The dashes are the same as the dots, except that the 
stationary solution has been normalized to match the maximum of /. For reference, the thin solid curve corresponds to initial conditions. 



5 ADAPTIVE REFINEMENT 

In the simulations with Gaussian and top hat initial conditions studied in previous section, the complexity of the phase-space 
distribution function augments with time, due to a well known effect of rolling up. In the top hat case, details are significant 
only in the external ring. This calls for local refinement, i.e. for increasing resolution only where needed: in principle this 
should reduce significantly the cost of the simulation. The gain obtained from local refinement in 2D phase-space might be 
of course quite questionable. For instance, such a gain is expected to be rather small in the Gaussian case, where complexity 
appears everywhere in the computing domain with approximately the same level of detail. However, refinement is expected to 
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Figure 16. The phase-space distribution function for top hat initial conditions, similarly as in Fig, 1121 
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Figure 17. Mean projected density and velocity dispersion, for the simulations with top hat initial conditions, similarly as in Fig, 1131 



be much more relevant in higher number of dimensions, especially in the cosmological case, where only a very small fraction 
of phase-space is occupied. 

Here, we examine a simple approach based on standard Adaptive Refinement Tree (ART) methods. Our goal is to only 
demonstrate that refinement is possible with the cloud method. This section is organized as follows. In S I5.1I we give a sketch 
of the method. Technical details are discussed in Appendix ^3] In § 15.21 we show how the method performs for the top hat 
initial conditions used in § 14.31 we start from low resolution initial set up with A x = A„ = 0.02, allow for two levels of 
refinement, and compare the results obtained to the high resolution simulations of § 14.31 with A x = A v = 0.005. 
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Fig. [TTI (continued - ). 



5.1 Method of refinement 

In the practical implementation that we test on the tree-code, the phase-space is decomposed hierarchically on a quad-tree 
at the moment of remap. Each cell of the quad-tree is associated to a cloud. If needed, the cell is splitted equally into four 
sub-cells corresponding to four sub-clouds. The process is performed as long as necessary, i.e. until the size of the corresponding 
(sub-)clouds obey some criteria based on local properties of the phase-space distribution function. 

At each successive level of refinement, we use Van-Citter algorithm to reconstruct the distribution function. We first start 
from the coarse level, where we reconstruct / the same way as for the fixed resolution code, as described in S 13.1.21 We then 
consider the set of clouds corresponding to first level of refinement as residues to reconstruct / more accurately by applying 
Van-Citter algorithm again on these clouds. And so on, until last level of refinement. 

We consider two criteria of refinement. The first one is based on the measurement of local curvature of the phase- 
space distribution function: local curvature is indeed a key quantity for preserving details in the distribution function: the 
cloud size should be always small compared to the local curvature radius. The second criterion is based on convergence of 
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Figure 18. Energy conservation for the simulations with top hat initial conditions, similarly as in Fig. 1141 

the reconstruction at successive refinement levels: clouds where / is poorly reconstructed are refined. We shall see that in 
practice, both curvature and convergence criteria give similar refinement structure. 

In our implementation of refinement, even when splitted into smaller sub-clouds, clouds at a given refinement level are 
thus kept since their sub-clouds counterparts are considered as residues. In principle, it would be advisable from the dynamical 
point of view to remove these clouds from the hierarchy, i.e. to set their mass to zero. Indeed, because of the structures in the 
gravitational potential brought by the small clouds, the quadratic approximation is expected at some point to become invalid 
for the larger clouds. We notice as well that one of our refinement criteria should rely as well on variations of the projected 
density, which account for deviations from quadraticity of the potential. Since our goal is just to demonstrate that refinement 
is possible with our method, we decided in the present work to put aside removal of coarser level clouds and refinement 
based on projected density. As it is, thus, our refinement procedure is quite improvable and will work only if the gravitational 
potential remains sufficiently quadratic at the coarse level, which should be hopefully the case if the number of refinement 
levels is not too large. 



5.2 Example: top hat initial conditions 

To check how our refinement procedure performs, we realized again simulations with the same top hat initial conditions as 
in § 14.31 but starting from low resolution initial conditions with A g = 0.02 and allowing refinement until high resolution is 
reached, A g = 0.005. This represents a significant increase in mass resolution, a factor 16, corresponding in total to a coarse 
level plus two levels of refinement. We performed two simulations, based on local curvature and local convergence criteria, 
respectively. All the parameters defining the simulations are the same as in §0] except for those depending on refinement, which 
are given in Appendices IB II and IB2I There is also a difference in the choice of 5 m i n (see § 13. 1.21 that we set to 5 m i n = 0.0005 
for simulations based on adaptive refinement with local convergence criterion. 

Figures 1191 1201 and 1211 show the phase-space distribution function at various times, t — 10, 40 and 100. The refined 
simulations seem to reproduce rather well the results of the full resolution simulation, up to t = 40. For t — 100, the refined 
simulations clearly differ from the full resolution one, even though they present the same level of detail. The two lower right 
panels on each of the figures show the refinement levels. As expected, the result obtained from local curvature criterion is 
very similar to that obtained from local convergence criterion. For reference, the upper right panel shows an estimate of the 
local curvature, more exactly the quantity C = max(|Ai|, | A2 1 ) intervening in Eq. 1B41 where Ai and A2 are the eigenvalues of 
the Hessian of the phase-space distribution function. Note that C is quite noisy despite our sophisticated weighting scheme to 
compute / and its derivatives, as it still captures some small scales defects, but it seems to be determined accurately enough 
to set up correctly refinement based on local curvature. 

Figures 12*2*1 and l2*3l show the projected density and the velocity dispersion for the snapshots t — 40 and t — 100. Their 
careful examination confirms the visual impression from Figs. 1201 and 1211 the refined simulations reproduce well the fine 
features of the full resolution one for t = 40, but not for t = 100. 

Finally, Fig. 1241 examines energy conservation for various simulations. Refined simulations present the same behavior, 
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whether local curvature or local convergence criterion is used: excellent energy conservation up to t ~ 40, then a nearly 
linear increase of energy with time with a significant but still reasonable final violation of energy conservation of the order of 2 
percent. The curves corresponding to the high and low resolution simulations of § 14.31 performed with Lucy deconvolution, are 
displayed for reference, but it would be more fair to make the comparison with fixed resolutions simulations performed with 
Van-Citter. Indeed, in our current implementation, simulations using Van-Citter deconvolution conserve energy less well than 
those using Lucy: on Fig. 1241 the high resolution simulation using Van-Citter presents a similar behavior to our simulations 
with refinement, except that the final violation of energy conservation is twice smaller. The low resolution simulation does 
significantly worse. First it conserves well energy, up to t ~ 10, then energy increases rather fast up to t ~ 40 where it reaches 
a plateau corresponding to a global violation of energy conservation of the order of 2 percent. 

This behavior is related to the fact that the positivity of the distribution function is not warranted in Van-Citter algorithm 
and that we perform a truncature of the tails of /: if / is smaller than 5 m i n , it is set to zero. If this truncature, which is 
meant to keep the computing domain finite and / as positive as possible, was not present, energy conservation would be by 
construction much better. The cut-off of the negative contributions of / implies that energy increases with time. The effect 
is all the stronger since resolution is low. Refinement is expected to improve energy conservation compared to low resolution, 
at least during some time, but since it keeps coarse and intermediary levels of refinement, it is not expected to do as well as 
the full resolution simulation, hence the result observed in Fig. 1241 

However, a two 2 percent violation of energy conservation is probably not enough to explain the significant disagreement 
at t — 100 between simulations with refinement and the full resolution one. There is indeed at least one other effect intervening, 
as discussed in last paragraph of § 15.11 namely that at the coarser levels, local quadraticity of the potential is likely to be 
violated during runtime. This can produce cumulative effects that show up only after a large number of orbits. 

Despite the limitations of the current implementation, we think having demonstrated here that the cloud method is 
compatible with adaptive refinement. 



6 PERSPECTIVES AND EXTENSION TO A HIGHER NUMBER OF DIMENSIONS 

We demonstrated the ability of the cloud method to solve the Vlasov-Poisson system in 1 dimension. An interesting feature 
of this method, is that its structure allows a simple generalization to the n-dimensional case. We will show in the next section 
that the n-dimensional cloud equations can be solved quite easily, as a natural generalization of the one dimensional equations. 
It is also quite obvious that the tree-code method that we used to solve the Poisson equation can be extended to a higher 
number of dimensions. 

6.1 n-dimensional equations. 

Provided that we keep in mind the results obtained in 1 dimension, the generalization of the cloud equation to n dimensions 
is straightforward. In the 1-dimensional case, it was shown that a general solution cloud be written as a function of second 
order polynomials in (xi...x n ,vi...v n ). We will adopt this result and show that a closed system of equations can be obtained 
for a n-dimensional cloud. To simplify the writing of the equations, we define the vector u such that: 

xt < i < n, 
Vi n < i < In. 

Using this notation, we are now able to write the generalized Vlasov equation in Lagrangian coordinates for a cloud of density 
/(u, t) evolving in the potential <f>{u, i): 

|-/(u, t) + Af( U) t ) u l+n n(i) - A ( U) t ) tt^-/(u, t) r,(i) = 0, (92) 

Ot OUi OUi OUi + n 



with the following definition for the function n: 

„(<) = / 1 ° <l ^ n ' (93) 

We will look for a general solution in the case of a quadratic potential. An extrapolation of the 1 dimensional case suggests 
that the general cloud equation can be written as a function of a quadratic polynomial. In this case: 

<t>(u,t) = OHj r)(i)r)(j) wuj, (94) 

/(u,i) = F(X ijUiUj ). (95) 

Note that the coefficients a^ and A^ depend on time, in general, ot,ij = 0Hj{i) and Xij = \ij(t). 
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Figure 19. Top hat simulations at t = 10. Left panels: phase-space distribution function for the full resolution simulation (top), the 
simulation with refinement based on local curvature (middle) and the simulation with refinement based on local convergence (bottom). 
Top right panel: the local curvature measured in the simulation with refinement based on local curvature, as explained in the text of 
§ 15.21 Two bottom right panels: levels of refinement in the simulations with local curvature criterion (middle) and local convergence 
criterion (bottom). Darker regions corresponds to higher level. 
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Figure 20. Same as in Fig. ED Dut f° r * = 40. 
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Figure 21. Same as in Fig. QUI but for t = 100. 
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Figure 22. Mean projected density (left panels) and velocity dispersion (right panels), as functions of position, for top hat simulations 
with refinement based on local curvature (top row) and on local convergence (bottom row), compared to full resolution (middle row), at 
t = 40. 
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Figure 23. Same as Fig. \22\ but for t = 100. 
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Figure 24. Energy conservation in top hat simulations with refinement, compared to fixed resolution simulations. The solid and dotted 
thin curves correspond to simulations with refinement based on local curvature and local convergence, respectively. The lower thick curve 
and nearly superposed dashed one correspond respectively to the high and low resolution simulations of § 14. .'{I These simulations use 
Lucy deconvolution. The upper thick and dashed curves correspond to similar simulations, but that were performed using Van-Citter 
reconstruction. 



By inserting this expression for /(u, t) in the Vlasov equation, Eq. 19211 . we will obtain a closed system of differential 
equations for the coefficients Xij(t). We will now illustrate the relevant calculations step by step. By taking each term in the 
Vlasov equation from left to right we have: 



First term in the Vlasov equation: 



Second term in the Vlasov equation: 



!'<»•*> 



dXi, 

dt 



UiUj. 



dui 



f(u,t)f](i) = 2\ijT)(i) uj 



Thus, 



Which can be rewritten: 



- — /(u, t) u i+n r)(i) = 2\j r/(i) UjU i+n . 



d 
— f(u,t) Ui+nVW = 2A (fc _ n)i rj(k-n) u 3 u k . 



(96) 
(97) 
(98) 
(99) 
(100) 



The last term of the Vlasov equation can be evaluated using the same method. We arrive at the following result: 

d d 

jr—<p(x,t)- /(u,t) T](i) =4<f> ik A( fe+ „)j 77(i)j?(fc) UjU k . 

OUi OUi+n 

We are now ready to tackle the Vlasov equation itself. We have to write that each coefficient of the second order polynomial 
is zero. For each second order term muj we have to consider the contribution from the pair of indices (i,j) and (j, i). In the 
case i = j, there is only one contributor, but we will still add the symmetric term, this will result in a factor of 2 on the left 
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side, and since the right side is equal to it does not change the equation. Thus by adding the contribution of each term in 
Eqs. 196L 1991 and IIUOL and taking into account the symmetric term, we find that: 

d\ 

-57 i + \i-™)3 V(i -n) + \j- n )i v(j - n) - 2(/> ki \ k+n)i T)(k)r)(i) - 20 fej A (fe+n)i r)(k)r)(j) - 0. (101) 

6.2 Applications of the method. 

Considering the current capabilities of the largest computers it seems clear that the 2 dimensional case could be undertaken 
rapidly, and that in astrophysics it would offer a new look on a few interesting problems, like for instance the dynamics of 
galactic disks. Due to the particular interest of the 2D case we present the detailed cloud equations in Appendix[U] Due to 
the high dimensionality of phase space in the 3 dimensional case, it seems that trying to integrate directly the Vlasov-Poisson 
system using the cloud method may be too costly. However, it is clear that a most interesting case, namely the cold dark 
matter model in cosmology, may not require the same amount of resources and could be solved using an appropriate version 
of the cloud method. It is important to notice that in this case the density in the 6 dimensional phase space is adequately 
represented by the extended folding of a 3 dimensional sheet having a nearly constant density. In particular, this case may 
not require the sophisticated Lucy deconvolution scheme and the general remapping technique. 
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APPENDIX A: SAMPLING OF THE DISTRIBUTION FUNCTION FROM A SET OF CLOUDS 

Al Arbitrary set of clouds to a regular grid 

Assume we want to compute fij — f(Xi,Vj) where (xi,Vj) are sampling points on a square grid of spacing A, for an ensemble 
of clouds with arbitrary axis ratio and orientation, using either Eq. I|59fl or Eq. I|66fl . Since the clouds are rather extended and 
overlap significantly, the calculation of fij, is expected to be rather CPU time consuming, in general, but simple. The natural 
way of performing it is simply to project each cloud independently on the grid. The only difficulty is to find the subsample of 
grid points where a given cloud contributes. To do that, one can estimate the projected size of the cloud on x and v axis: 

T-, _ Jtmax / P2 t-> Amax / PO / . -. \ 

"max,i — p y -^, -rtmax,u — — A/ — ; , l Ai j 
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where /3o and /?2 are given by Eqs. 1221 and 12 U . respectively. This defines the rectangle where the ellipse corresponding to the 
cloud is inscribed. To speed up the calculation, one can find, for fixed values of x £ [--R m ax,i,-fimax,i] in the cloud coordinate 
frame, the segment [vi, V2] given by the intersection between the ellipse covered by the cloud and the vertical line of abscissa 
x, that we do not feel necessary to write here for simplicity. 



A2 Round set of clouds to a grid: fast convolution algorithm 

When one considers a round set of clouds set up on a regular grid, it is possible to use a simple trick to speedup the calculation 
of fi t j compared to the simple algorithm described in § IA1I by noticing that the convolution of a function by a Gaussian can 
be factorized in that case: 

/ dx* dv* f(x + x*,v + v*) exp -- ( ^T + ^J ) = / dx * exp [^2^W) / dV /( x + x '*> ^ + "*) exp r2¥ ' ^ A2 ^ 

In practice, that means the following: assume for simplicity that the center, (xq,Vq), of each cloud, n, coincides exactly with 
a grid sampling point (x% n , Vj n ), and furthermore that cloud inter-spacing, A g , is a multiple of the target mesh inter-spacing 
A. With these assumptions, first create an array /**, which is zero everywhere, except at cloud positions: 

ft: )jn = M n /V, (A3) 

where M n is the mass of the cloud n, and V = Ag. Then perform an operation equivalent to the convolution in velocity space. 
For doing that, it is just necessary to propagate vertically the initial values f** j n using the weights given by function G(x, y), 
on the columns of the array f** which contains non zero elements: 

The last operation consists of convolving along x axis: 

kj = mm ^ fLjGfa - x%,o). (A5) 

Here we used the separability of function G(x,y): G(x,y) = G(x,0)G(0,y)/G(0,0). We made the approximation, to speed up 
the calculation, that the truncation is not x 2 + v 2 < -R^axi but rather \x\ < i? max and \v\ < -R ma x, but this should not have 
any significant consequence if the value of i? max is large enough, e.g. -R max = 4R as used in this paper. For simplifying the 
argumentation, we introduced the arrays /," and /*,■ with the same size as the target mesh size, but one can see that only 
/* is necessary, and it can be reduced to an array of dimension (N C:X ,N gtV ), where N c ,x and N gtV are the size of the cloud 
mesh along x axis and of the target mesh along v axis, respectively. Typically, the first step of the calculation, Eq. (IA4II . takes 
2Aftotiimax/A operations, where JVtot is the total number of clouds, while the second step, Eq. 1A50 . takes 2A / a ites-Rm a x/A, 
where JVsitcs is the total number of target grid sites. The total number of operations, (JVtot + A r s itcs)2_R ma x/A, is to be compared 
to 7rJV s itcs(i? ma x/A) , as expected from the method presented in S IA1I a tremendous gain in time. Indeed, if we use Eqs. IJA40 
and 1A5II to compute ft in Lucy or Van-Citter algorithms, we see that we gain a factor of order of (4/7r)(i? max /A) 2 ~ 10 for 
the deconvolution and our parameters choice (J? max = 47?, R = A/\/2). As a result, we used this fast convolution method in 
the CM code to speed up reconstruction and the method described in § lAll in other cases, when the clouds have arbitrary 
shape and orientation. 

In principle, we could have used this efficient implementation as well for our tree-code. However, this program is designed 
for adaptive mesh refinement (§|5J. It would be rather involved, algorithmically, although possible, to adapt this method to an 
unstructured grid such as the one used in our refinement procedure. Instead, for all the simulations presented in this paper, 
which are done with the tree-code, we used the very general but much slower method described below, at a significant cost of 
CPU time. 

Note finally that Eqs. <A4t and IA5t apply to Eq. 15911 . but the method can easily generalized to Eq. 1661 . 



A3 Arbitrary set of clouds to arbitrary point: the quad-tree algorithm 

In general, if one aims to estimate f(x, v) at an arbitrary point of phase-space, it is necessary to find rapidly the clouds 
contributing to this point. There is a very general and standard method to do that, based on hierarchical decomposition of 
phase-space on a quad-tree structure, until there is zero or one cloud center per cell of this tree. The list of clouds contributing 
to point (x,v) is constructed by walking into this tree, from root to leaves. To build the tree properly, one has take into 
account of the extension of the clouds. While the tree is constructed, for each cell it contains, defined by some coordinate 
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range, [(xi, Vi), (X2, V2)], one computes the effective potential range C e s = [(s e ff,i, v e ff,i), {x e s,2,v B s,2)\, where 

Xeff.i = min (xq - Rln^x) , (A6) 

clouds i in the cell V / 

x c g,2 - max ( Xq + i4ax,! J , (A7) 

clouds i in the cell \ / 

VeB.i = , , min (vjj - i4 K| J , (A8) 

clouds 1 in the cell V / 



«eff,2 = max I «g + i4»,« ) • (A9) 

clouds i in the cell \ ' 

In these equations, J?J„ axl and J4 ali „ are given by Eqs. lAH for each cloud belonging to the cell, {xq,Vq) € [(xi, Vi), (X2, v^)]. 
To construct the list of clouds contributing to point (a;, v), one starts from the root and walks down into the tree. During the 
walk, only cells verifying (a;, v) £ C e g are opened, i.e. decomposed in four sub-cells, until convergence is achieved, i.e. when 
the cell contains zero or one cloud. This algorithm is quite fast: the loss of speed compared to the trivial method explained 
in S IAll is of the order of IniVtot, where iVtot is the total number of clouds. It has the advantage of being very general since it 
allows to compute f(x, v) at any point. This explains why we adopted it from the very beginning when we started to develop 
our tree-code, in order to keep this latter as flexible as possible. 



APPENDIX B: REFINEMENT PROCEDURE: TECHNICAL DETAILS 

This appendix deals with technical issues of our refinement procedure. In S IB1I we explain in detail how the refinement structure 
is implemented. In particular, as in standard ART methods, to keep the method stable, we make sure that transitions between 
refinement levels are not too abrupt. We also take into account the large extension of the clouds. In § IB2I we describe the 
way we set our criteria of refinement, namely the one based on local curvature and the other one based on local convergence. 

Bl Principle of refinement 

Let L be the level of refinement, equal to L = for the coarse level. For simplicity, we assume that the coarse level is a grid 
with fixed inter-spacing, A^ = A v — A g . In order to avoid instabilities, function L(x, v) should vary smoothly in phase-space: 
abrupt transitions between two levels |Li — L^\ > 1 are forbidden. Also, special care has to be taken of the large extension of 
the clouds. For example, assume that we split a cloud of size R = %/2/2A g with cut-off radius 7? ma x = 47? into 4 twice smaller 
clouds. Since the cloud influences remote regions beyond the cut-off scale of the 4 sub-clouds, refinement should be performed 
as well in these regions. Fig. IB1I shows how to handle this problem and to naturally preserve smoothness of function L(x,v). 
While following the rules dictated by the procedure detailed in Fig. IB1I we use a criterion based on local convergence or/and 
on measurement of local curvature to decide if a cloud has to be refined or not, as explained later in § IB2I 

Once a full hierarchy of clouds and sub-clouds has been created, we proceed as follows to reconstruct the phase-space 
distribution function: 

(o) Setup coarse level, L — 0: compute ff from the old set of clouds (prior to remap) at sampling points (xq , Vq ). Then 
use Van-Citter algorithm to compute the masses, Mf, of the coarse level clouds, as explained in § 13.1.21 

(i) If maximum level of refinement is not reached, increase L by one, otherwise stop. 

(ii) Compute ft at the positions of clouds of level L, (xq ,Vq ) , from the old set of clouds. Now, consider clouds of level 
L like residues and use Van-Citter algorithm to compute the masses of these clouds, M/' . To estimate ft i n Eq. 1631 , clouds 
of lower levels £ < L — 1 contribute, but with fixed masses computed previously. Start again with step (i). 

In step (ii) , once the cloud masses have been determined up to some level, they remain unchanged while determining the masses 
of clouds of upper levels. It is indeed expected that the corrections brought by residues improve the reconstruction of the 
distribution function within the corresponding refinement region, without contaminating the remaining part of phase-space. 

In our implementation of refinement, we use Van Citter deconvolution algorithm. Alternatively, one might use Lucy 
method. In that case, clouds with L > 1 cannot be considered anymore as residues, and, due to the multiplicative nature of 
the Lucy algorithm, we would have to reiterate at all levels to compute the overall distribution of masses, Mf, £ < L. This 
would complicate considerably the algorithm, with no warranty of easy convergence, so we decided for this work not to test 
such refinement scheme. 

To compute fi in Eqs. 163H and 1641 from the new set of round clouds, one can still use Eq. I59H . as discussed in end of 
§ 13.1.31 It now simply reads 

f(x,v) = JT ±- t ]T Ml G t (x ~ x%',v- v^), (Bl) 

1=0 i 
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Figure Bl. Sketch of our refinement procedure, for our choice of cloud size, R^, compared to local cloud inter-spacing, Ag/2^: 
2 Rl I A g = v2/2, which implies a 4i?£ cut-off at 2 i i? max /A g = 2\/2. Suppose we decide to refine the right most cloud of upper panel, 
of level L — 1. This cloud is supplemented with four twice smaller sub-clouds of level L, labeled by 3. It however influences a region 
of radius -R max /2 i— x , represented by the black arc of circle. To enforce our refinement criterion on all the region influenced by this 
cloud, we must add supplementary layers of sub-clouds around it, labeled by 2 and by f . We can stop adding corrections at sub-clouds 
labeled by one, since their extension covers the zone of influence (light arcs of circle). For easy implementation of regularization IB2I . 
we create as well ghosts clouds of zero mass designed by the "G" letter in order to compute correct weights for the refined level. With 
this refinement procedure, we see that the cloud to be refined has to be sufficiently far away from ghosts. The closest possible ghost for 
the rightmost cloud at level L — 1 is the left most one, which was labeled by a G surrounded by a circle. 
Our refinement can be performed recursively to higher level, as illustrated by lower panel, by obeying the following simple rules: 

(1) it is forbidden to refine ghosts; 

(2) clouds labeled by I can only be refined in ghosts; 

(3) clouds labeled by 2 can only be either refined in ghosts, in sub-clouds labeled by 1 or in sub-clouds labeled by 2 according to 
their distance to closest ghosts, as illustrated by the upper and lower panel; 

(4) clouds labeled by 3 can be fully refined. 

Note that, by construction, the region labeled with 3's can be refined to arbitrary level Loo, while keeping smooth transition between 
coarse level and Loo- 

where V — V°4~ e = A :c A u 4~ f and function Ge is the same as function G but with a radius Ri = 2~ R. To estimate //" 
from the ancient set of clouds, one can similarly adapt Eq. 1166(1 . which now becomes 



f{x,v) 



2.^ \/i 
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)Mi 



V 1 



Ei Gl(2 



(B2) 



<-G >' 



We see in Eq. HB2II that the interpolation is performed at each level of refinement separately. In order to have proper 
normalization, we need to add at the border of each refinement level a layer of ghost clouds with zero masses but contributing 
to the weight in the denominator, as described in Fig. IB1I 



B2 Criteria of refinement 



Our refinement is by essence set up to follow details of the distribution function where needed. Consider the distribution 
function as a surface of equation z = f(x, v). The local curvature of this surface determines to which extent these details can 
be reproduced by our clouds of finite size. If the size R of the clouds is not small enough compared to the local minimum 
curvature radius 1Z C , we expect significant loss of details as well as aliasing effects. A first natural refinement criterion thus 
relies on the value of R/1Z C - However, it is expectable, equivalently, that deconvolution will have trouble to converge when 
R/TZc > 1, since details cannot be adjusted correctly by the clouds. This motivates an alternate criterion of refinement, based 
on local convergence of the reconstruction. The convergence criterion has also the advantage of being a quite natural extension 
of our approach: it guarantees by definition the quality of the reconstruction, at variance with the local curvature criterion. 
This latter indeed relies on the measurement of second derivatives of /, which can be quite noisy. However, as illustrated by 
§ 15.21 both refinement criteria give very similar results. 
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While constructing the hierarchy of clouds and sub-clouds, we proceed level by level. To create clouds of level L, clouds 
of level L — 1 are tested as follows: 

(A) Local convergence criterion: in that case, it is necessary to reconstruct the distribution function simultaneously with 
the refinement structure, using steps (o), (i) and (ii) described in § IB1I Simply, step (i) has to supplemented with the criterion 
on local convergence to create the new level of refinement, L: one sets a list of clouds of level L — 1 for which the reconstruction 
scheme did not converge, S^ _1 > S^ 1 in Eq. I63L after some fixed number of iterations, say 10 according to discussion in 
§ 13.1.21 In principle, the optimal calculation of S^ -1 as a function of L depends on the noise properties of the reconstruction, 
but it would go beyond the scope of the paper to analyze them in detail. We therefore consider a very simple way of setting 
the convergence threshold, 

Sc = S° c , L > 1, (B3) 

which ignores possible propagation of errors in the reconstructions from level to level, at the risk of accumulating small scale 
artifacts. For this reason, the value of 6° should in practice be slightly larger than discussed in § 13.1.21 to avoid unnecessary 
refinement due to fluctuations of the small scale noise. A good practical value is <$J? ~ 0.001. 

(B) Local curvature criterion: the clouds of level L — 1 for which the local curvature of the phase-space distribution function 
is larger than some threshold are refined. The local curvature is estimated at cloud centers from the old set of clouds (prior to 
remap), using the Hessian of the function f(x,v) given by Eq. 1B21 (we do not write it here for simplicity). Our refinement 
criterion is 

A \ 2 
-^ I max(|Ai|, |A 2 |) > F curvaturc , (B4) 

where Ai and A2 are the eigenvalues of the Hessian. A good practical choice of the refinement parameter is F curV aturc — 0.1. 

(G) Keep stability of the refinement: at positions of clouds of level L — 1 to be refined, check locally what was the refinement 
level L id in the old set of clouds and impose that L < L ia + 1. To do that, we consider the square of size A g /2 L_1 centered 
on (xq ~ ,Vq _1 ) and find with standard quad tree search the set of old clouds intersecting with it. L id corresponds to the 
maximum refinement level found in this set. For stability purpose, we should have L < L id + 1. During runtime, we indeed 
do not expect the creation of more than one level of refinement per remap. If that happens, it must be due to some spurious 
small scale artifact created at the moment of previous remap, that we do not want to propagate furthermore. 

There is a subtlety in our algorithm that we have to mention now. Indeed, the criteria explained above are applied to a 
discrete set of sampling points. To make sure that this sampling catches all the features of the distribution function, one has 
to examine locally conditions (A) and (B) at all levels up to I/ id, by creating the corresponding sub-clouds. Due to the large 
extension of the clouds and the way L id is computed, a large number of sub-clouds are in fact unnecessary. At the end, the 
hierarchy of clouds has to be "cleaned up", to keep only the relevant refinement regions, i.e. regions verifying criteria (A) 
or/and (B). After being cleaned up, the new set of clouds and sub-clouds is ready for successive deconvolutions at various 
levels, which are thus performed twice when convergence criterion is used: once to construct the hierarchy of clouds and once 
again after this hierarchy has been cleared out from unnecessary sub-clouds. 



APPENDIX C: EXPLICIT 2 DIMENSIONAL CLOUD EQUATIONS 

'' /Al I_4 11 A 31 -4 21 A 4 i=O, (CI) 

4</>12A32 — 4 022A42 = 0, (C2) 

+ 2Ai 3 =0, (C3) 

ilf +2A 24 =0, (C4) 

—- 2 0i2 A31 — 2 022A41 — 2 0n A32 — 2 02i A42 = 0, (C5) 

^-+Aii-2 0nA33-2 02iA 4 3=O, (C6) 



dt 


— 4<Pli 


0IX22 
dt 


-4 012 




dA 33 




dt 




a A 44 



46 C. Alard & S. Colombi 



^ + A12 - 2 0i 2 A 33 - 2 22 A 43 = 0, (C7) 

dA4i 



+ A21-2^1lA 3 4-2^2lA 4 4=0, (C8) 

+ A 22 -2< ? ii 2 A34-2 022A 4 4 =0, (C9) 

+ A 23 + A14 = 0. (CIO) 



dt 

dt 

d\ i3 
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